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SEPARATION OF LAPLACE’S EQUATION* 


BY 
N. LEVINSON, B. BOGERT anv R. M. REDHEFFER 
Massachusetts Institute of Technology 


ABSTRACT 


The following results are established in this paper: 

(1)** For the Laplace equation A@é = O in curvilinear co-ordinates (u, v, w) in Euelid- 
ean space to be directly separableT into two equations, one for S and one for Z, when the 
solution is @ = R(u, v, w)S(u, v)Z(w) with fixed R, it is necessary and sufficient that the 
surfaces w = constant (1) be orthogonal to the surfaces wu = constant, v = constant and 


(2) be parallel planes, planes with a common axis, concentric spheres, spheres tangent at a 


common point, or one of the two sets of spheres generated by the co-ordinate lines when 


bicireular co-ordinates are rotated about the line joining the poles or about its per- 
pendicular bisector. 

(II) We have R | always and only in the first three cases, namely, when the sur- 
faces w = constant are parallel planes, planes with a common axis, or concentric spheres. 

(111) In these three cases, but only these, the wave equation separates in the sense 
RSZ, and hence, tor the wave equation, R = 1 automatically. 

(IV) For further separation of the equation found above for S, when S = X(u)Y(v) 
so that the solution is now RN YZ, it is necessary and sufficient that the co-ordinates 
be toroidal, or such that the wave equation so separates, or any inversions of these. 

(V) The co-ordinates where the wave equation so separates, that is, admits solutions 
RX (uw) Y(v)Z(e), are only the well-known cases where this happens with Rk = 1, namely, 
degenerate ellipsoidal or paraboloidal co-ordinates (but see See. 8.2). 

(VI). In these cases, but only these, R = 1 for the Laplace equation too. 

VII) Co-ordinates for RSZ or RN YZ separation of the Laplace equation have the 
group property under inversion. 

VIII) In all eases R can be found by inspection of the linear clement. 


INTRODUCTION 


0.1. The problem. Separation of variables is one of the simplest and most frequently 
used methods of solving partial differential equations subject to given boundary condi- 
tions. In this method, with which we assume the reader to be familiar, the surface over 
which the boundary values are specified must coincide with one of the co-ordinate 


*Received June 28, 1948. The work of Redheffer in this paper was supported in part by the Signal 
Corps, the Air Materiel Command, and the O. N. R. 
**The se parate re sults are numbered for ease of reference. 


See See, 0.2 
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surfaces, and hence any restriction on the co-ordinates which can be used is also a re- 


striction on the physical situations to which the method can be applied. It is natural 


il 
to inquire how we shall « 
can be solved by separation of variables. Such is the type of problem with which the 


letermine all co-ordinate systems in which a given equation 


present paper is concerned. In a later paper we give a brief history of the subject and 
a@ more general discussion of separation methods. 

Partly because it includes the wave equation in problems of this type (ef. Sec. 4.3), 
and partly because the results are still unknown, the emphasis here is on separation of 
variables in Laplace’s equation. We confine the discussion to three dimensions as the 
case of greatest practical interest. Also this loss of generality, which is believed to be 
inconsequential, enables us to put the derivation in a form readily followed by the 
general reader whose interests are not primarily mathematical. Such considerations have 
weight because the problems treated are of interest to physicists as well as mathe- 
maticians. 

0.2. The manner of separation. Turning now to the question of the mode of separa- 
tion, we observe that a more general form than the XYZ generally considered will 
achieve the objectives desired. Specifically we may assume the form RXYZ, where R 
is a single fixed function’ of the three co-ordinates: R = R(u, v, w). If X, Y and Z form 
a complete set * and R is known, the operations can be carried out with no increase in 
complexity. It might be thought that.a partial differential equation would have to be 
solved to get R, but we shall see that this is not the case. 

Assuming the form RX YZ rather than XYZ, we have a slightly milder restriction 
on the solutions and we might expect, therefore, to find a larger class of permissible co- 
ordinates. Such is indeed the case; with toroidal co-ordinates, for example, the equation 
is known to separate in the extended sense but not in the restricted sense [11], and the 
same is true for the so-called Dupin cyclides [7]. 

A discussion of separation must consider the method of obtaining the solution, as 
well as its form. In most circumstances the solution is obtained by separating the 
equation in the literal sense; that is, we multiply through by a fixed function of the co- 
ordinates to obtain two sets of terms, one set involving w only, for example, and the 
other involving u and v only. The same procedure is then used for the u, v terms. In 
the present work we assume not only that the solution has the form XYZR, but also 
that the equation can be separated in this way. Such an assumption is restrictive; in 
particular we do not obtain the Dupin cyclides even though the solutions have the 
prescribed form. One of the categories of co-ordinates in [3] is also omitted. In a sub- 
sequent paper we consider the general case, with its relation to this and other special 


modes of obtaining the solutions. 
I. Sotutions R(u, v, w)S(u, v)Z(w 


SEPARATION OF THE EQUATION 


1.0. Laplace’s equation, and the second derivative terms. Consider the Laplace 
Pp quack 


operator under a transformation of co-ordinates from (x, y, z) in Euclidean three- 





1Feshbach has raised the question of co-ordinate systems for which the wave equation separates in 


this sense. 
2This means that a sufficiently well-behaved, but otherwise arbitrary, function f(z, y) can be ex- 


panded in the form 2 AqsXaY; ; similarly for other pairs (X, Z), (Y, Z). 
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dimensional space to curvilinear co-ordinates (u, v, w). Since the appearance of cross 
derivatives involving w makes separability impossible, we consider at first the case of 
orthogonal co-ordinates. It will be found later that the general case of oblique co- 
ordinates can be easily reduced to this one (Sec. 4.1). In the curvilinear co-ordinates 


(u, v, w), then, let the linear element be 
ds* = f? du’ + g° dv + h* dw” (1.0) 
where f, g, and h are functions of u, v, and w. Then [11] 


1 0 {gh ) ri) (Af ) 0 (f2 )] 
g = —| — (= — {=< —( 
. fgh E (4 bu) + dv\g 8) + dw \h O 


which becomes, if @ = R(u, v, w)S(u, v)Z(w), 


17D I Suu I S.. ie «i 
A = szn| } 9 + 7 8 a 
(1.1) 
. a * Z’ us 
+7 46,545,54746,4 B 
where 
7 _~ 2h, 1 (2) 


and similarly for G, and H, . Also 


F — I Ru + / 1 (2%) R., (1 3) 
- Fes foh\f/,R ; 
and similarly for G, and H,. 
For separability of A*@ = 0 in the form (1.1) we require that the function Z separate 
off into an ordinary differential equation. That is, there must exist a function A*(u, v, w) 
such that when Aé@ is multiplied by A*, only the variable w appears in the coefficient 
of Z"/Z and Z’'/Z, and w does not appear in the coefficient of any other differentiated 
terms. Moreover A°(F, + G, + H,) must break up into the sum of a function of w 
only and a function of u, v only. 
Using the condition on the coefficients of the terms involving second derivatives of 
S and Z, and bearing in mind the definition (1.0) of f, g and h, we see that the linear 
element must have the form 
ds? = A*[F’ du’ + G’ dv? + H’ dw’, (1.4) 
where F, G are functions of u, v only and H is a function of w only. 
1.1. Permissible changes of variable. If we replace w by a function of w, and the 
u, v co-ordinates by new co-ordinates which still do not involve w, we may put (1.4) 
in the form 
ds’ = A*[F? (du* + dv’) + dw’] (1.5) 
where F is a function of u and v alone. Here we have written A for the new value of A 


and u, v, w for the new values of u, v and w. In addition, we have used the fact that any 
element F? du? + G’ dv’ can be mapped conformally on a plane. 
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The w transformation is certainly permissible, as it amounts merely to a change of 
seale for the w co-ordinate. Hence it does not alter the geometrical configuration of the 
co-ordinate surfaces. When we say that a change of variable is permissible, in this con- 
nection, we mean that the geometric properties relevant to separation of Laplace’s 
equation are essentially unchanged. Hence proof of the existence of a certain geometric 
property in the transformed system shows that the same property was present in the 
original system. 

That the change of u, v co-ordinates is permissible in this sense follows from the 
fact that the new linear element will have the form (1.4). Hence the condition that 
the equation should separate is still satisfied (note that Eq. (1.8) will also persist in 
form). Moreover, since every linear element of the form (1.4) corresponds to a triply 
orthogonal system, we see that the equations used later, (2.0) and (2.1), will continue 
to be valid in the new system. From these and similar remarks concerning changes of 
co-ordinates we conclude finally that it suffices to specify the surfaces w = constant. If 
the equation separates in our sense for a particular system of u and v surfaces orthogonal 
to these it will separate for every other such system. Conversely, if a change of u and 
v co-ordinates leads to a certain set of u and v surfaces, and it is then found that the 
w surfaces must have certain special properties, these properties will in fact persist for 
every choice of the u or v surfaces. To simplify the analysis these two operations, changing 
the w scale and changing the wu, v co-ordinates, will be repeatedly used in the ensuing 
discussion. 

1.2. The first derivative terms. Turning now to the coefficients of terms involving 
first derivatives, we use the fact that AH, must be a function of w only to find, by virtue 
of (1.2) and (1.5), that 

ul 2 kt + z -(F°A), | = h,(w), 
A’ R d ‘ 


which simplifies to 
2k A 
= + = fi,(w). 
R A 
We integrate with respect to w, noting that the constant of integration may be an 
arbitrary function of u and v, and we take the exponential of each side, to find finally 


RA = h.(w)F3(u, v). 1.6) 
Since it is permitted that R involve u, v, and w we may always modify F in such a way 


that 
BA = 4%. Lz) 


Thus, the term /, in (1.6) may be absorbed in Z and F; in S. Equation (1.7) leads to 
the result: on¢ may determine R explicitly by putting the linear element in the form (1.5), 


and then taking R = | A’ 

It may be verified that the coefficients of S,/S and S 
u and v only, without any new condition (actually they are zero). From separation of 
7), and an extra condition, 


S are already functions of 


the equation, therefore, we have only (1.5), (1. 
A*(F, + G, + Hz.) = hy(w) + Flu, vr). (1.8) 


It will be found that (1.8) is a consequence of the Euclidean character of the space. 








Se 
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USE OF PROPERTIES OF EUCLIDEAN SPACE 
2.0. The functional form of A. If the space is Euclidean then the linear element 
(1.0) satisfies the relations [2], [6], 


| Gufe _ hfe _ g (2.0) 
g h 


B 4 || + JG = 2. (2.1) 
fi” ted h 


We also have those relations obtained by simultaneous cyclic permutation of (fgh) and 
(uvw). These equations, which state that the Riemannian curvature of the space is zero, 
give a necessary and sufficient condition that the co-ordinate system be imbedded in 
Kuclidean space. 


In terms of (1.5), the relation (2.0) containing g,,, becomes 


foe me go (2.2) 


after simplification. Integrating with respect to w, as in the derivation of (1.6), we find 
that (1/A), is a function of u and v alone. The same is true of (1/4), . Thus it follows 
that 1/A = t(u, v) + H(w), or 


A = 1/|t(u, v) + H(w)). (2.3) 


The two relations (2.0) used in the derivation are now satisfied, and in terms of the 
new linear element the third gives 


t., = F,t,/F + F.t,/F. (2.4) 


‘Turning now to the relations (2.1), we find that (2.1) as it stands leads, without 
detailed calculation, to an equation of the form* 
H’ + ¢,H* + 2c.H +c, = 0, (2.5) 
where the coefficients ¢ are functions of u and v only. Differentiating with respect to w 
we find that 
H'(H” + oH + e) = 0. (2.6) 
Since the solution depends on w only, while the coefficients depend on u and v only, we 
may give to u and v any constant values to get a linear differential equation for H with 
constant coefficients. This equation may be solved to give 
H = 0, wv, w’, e“, sinh w, cosh w, sin w (2.7) 
as the only possible values of H that are essentially distinct. Here the case H’ = 0 has 
been reduced to H 0 by absorbing the constant value of H in the function ¢. If H 
is not zero we substitute back into (2.6) and find that the coefficients, which we already 
knew were independent of w, are actually constant. This result will be needed later. 
In summary, we have found that the linear element must be of the form 
F?(du? + dv’) + dw’ 
= . — (2.8) 


ds i + Hy? imi 


‘The equation is written explicitly in (2.9). 
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with ¢ and F functions of wu and v only, while H is one of the functions (2.7). The only 


additional conditions are (1.8), (2.1), and (2.4). 
To use the relations (2.1) we substitute for f, g and h their values as given by (2.8). 


The equation (2.1) as it stands gives 


2 2 2 72 a fu 
fet fe REET EH. (B) + (2) (2.9) 
t+H (¢ + H) F /, F/, 

while the sum of the two others found by permutation leads to 

ti. +t, + F°H” 12 
he + ty = 22> ——— — 2F°H" (2.10 
J+ t, t+ H F°H 10) 
and their difference gives 
F, ae, 

— t,, = 24,5 — 24,3. 2.11 
bn = tne lu F lb. 5 (2.11) 


The system (2.9)-(2.11), slightly simpler than (2.1), is equivalent to it. These three 
equations and (2.4) give the necessary and sufficient condition that the metric (2.8) be 
for orthogonal co-ordinates in Euclidean space. The fact that ds* has the form (2.8), 
plus the relation (1.8), on the other hand, gives the necessary and sufficient condition 
that this co-ordinate system be one in which Laplace’s equation separates. Systems 
satisfying both sets of conditions, and only those, are the ones we are seeking. 

2.1. The second fundamental form. It is known from differential geometry [2], [6] 
that a triply orthogonal system is completely determined, except for its orientation in 
space, by the linear element ds’. Hence the second fundamental form (as well as the 
first) for each of our co-ordinate surfaces is specified by f, g and h of Eq. (1.0) and we 
might expect to find specific relations from which they could be computed. Such rela- 
tions actually exist [2], [7]; on the surface w = constant we have 


—dz-dN = ldw + mdudv+ndw’, 


where 
l= —ff,/h, m = 0, nN = —G9Go/h; (2.12) 
the results in other cases are obtained by cyclic permutation of u, v, w and f, g, h. 
In particular for the linear element (2.8) we obtain that 
—- F°H’ ; 
—dzx-dN = — 3 (dué lv’). (2.13 
sail “tn ©? (2.18) 
On the other hand setting dw” = 0 in (2.8) gives 
dz-dz = —— Fr ; (du? + dv’) (2.14) 
(t+ H) 


for the first fundamental form. Since the two forms (2.13) and (2.14) are proportional, 
we know that the surfaces w = constant must be spheres or planes [13]. Computing 
the Gaussian curvature as the ratio of the two discriminants d’/D* (Ref. [13]) we find 


the radius of the sphere corresponding to a given value of w: 
radius = 1/H’(w). (2.15) 


This result will frequently be used in the ensuing investigation. 
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SEPARATE EXAMINATION OF CASES 

3.0. The case // = 0. For each value of H in (2.7) we obtain a new set of relations 
for F. These equations are not easy to solve as they stand, and our procedure will be 
to seek a change of u, v variables that will reduce them to simpler form. By the dis- 
cussion of Sec. 1.1 we know that any restrictions on the w surfaces obtained after the 
change must have been valid before it as well. 

In case H = 0, as assumed here, we see by (2.13) that 1 = m = n = O and hence 
the surfaces w = constant must be planes [13]. On one of these planes let us pick the 
u, v co-ordinates so that the linear element takes the Cartesian form du? + dv’. This is 
possible, since the surface is a plane. By comparison with (2.8) when dw* = 0 we see 
that the present u and v co-ordinates make F/t = 1. Also since both ¢ and F are inde- 
pendent of w, we now have F/t = 1 for all values of w, not merely for the constant 
value first selected. 

With this procedure Eqs. (2.10) and (2.11) become 


9° 


am + tas _ 2(t, + t.)/t, 


tun — tye = A(t; — t)/t, 
whence, after adding and dividing by ¢, , 
buu/ty = 2t,/t 
with a similar result for ¢,, . Integrating twice we find that 
1/t = uf,(v) + fr(v), 


and a corresponding result. with wu and v interchanged. The two together show that 1/¢ 
must have the form a + bu + cv + dw, with a, b, c, d constant; Eq. (2.4) tells us that 
d = 0. The remaining conditions (2.9) and (1.8) are now satisfied, so that there is no 
other restriction. 

If b = c = 0 the surfaces w = constant represent parallel planes, as we see by com- 
parison with the linear element for Euclidean co-ordinates. If b or c is not zero, however, 
we have ds? = du? + dv? + u? du’ or, renaming the variables, ds? = dr? + dz’? + 1° dé’. 
This shows that the w = constant surfaces are planes with a common axis. 

3.1. The case 7 = w. Turning now to the case H = w we find from (2.10) that 


(tu + te )(t + w) = 2 + & + F’). 


Since this is an identity in w the coefficient of w must vanish. Consequently ¢,, + t,. = 0, 
and hence also #2 + é + F’ = 0. This is impossible since it makes F = 0. 
3.2. The case H = w’. Next if H = w’ we find that 


(F./F), + (F./F), = 0, (3.1) 
ti. + t,, = 4F’, (3.2) 
tu +t, = (4 + t)/t, (3.3) 


by equating to zero the coefficients of w‘, w’ and 1 in (2.9). Nothing new is obtained from 
(2.10), and hence these with (2.4) and (2.11) are the only conditions. Equations (3.1) 
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and (3.3) are equivalent to the statements that log F and log ¢ are harmonic functions. 
By simply writing the Gauss equation for curvature in terms of the first funda- 


mental form, we find that for any F satisfying (3.1), this curvature vanishes for the 


surfaces with ds* = F(du*> + dv*). These surfaces therefore are developable, and we may 
introduce a new set of u, v co-ordinates which will make F = 1. When this is the case, 
Eq. (2.4) tells us that ¢,, = 0, so that ¢ = a(u) + B(v). Equation (3.2) now reduces to 
a’ + B” = 4, which gives a” = 2+ ¢, B” = 2 — c withe a constant. It follows that 


a =u, =v’, after making a change of variable, if necessary, to eliminate the arbitrary 


constants. Our linear element now takes the form 
9 due + di 4 + dw” - 
ds = . . > ( 4) 
(uw +v+ wu’) 


after ¢ has been given its value as determined above. We observe from (7.1) and (7.2) 


that the linear element (3.4) is the one which would be obtained by inversion of Car- 
tesian co-ordinates in the origin. Because the linear element is sufficient to determine 
the co-ordinate system completely, as noted in Sec. 2.1, it follows that the co-ordinate 
surfaces, as well as the linear element, will be the same as those which would be obtained 
by the inversion described. In particular the surfaces w = constant must consist of a 
plane and a set of spheres all tangent to it at one point. This result, incidentally, can be 
obtained directly from (3.4), as in the examples considered below. 


3.3. The case H = e”. Turning now to the case // = e” we substitute in (2.10) and 
equate to zero the coefficients of 1 and e”, respectively, to find: 
t(tu. + te.) = At. + &), 
liu te —2k't 
which gives 
t+tie+Ft =0. 3.5) 
Equation (3.5) implies that ¢ is zero. We have then from (2.8) that 
: F*(du- + dv’) dw . 
ds = : = + 3.0) 
The other relations are now all satisfied, if F is suitably restricted. 
By the second fundamental form the radius of the sphere w = constant is e “. The 
distance from w = & to the sphere along u = a, v = bis also e “, by Eq. (3.6). Thus 
the spheres are concentric with w = © as the common center. 


3.4. The case H = sinh w. When H = sinh w, we substitute into (2.10) as usual, 
then put everything in terms of sinh w by using cosh” w = 1 + sinh’ w, and finally 


equate to zero the coefficients of 1 and sinh w. We thereby obtain 


“(t,. + t..) = 2Ai+te+F), (3.7) 
i, +t, = —Hs. (3.8) 

Together these relations show that 
P4647 +4 Ff = @ (3.9) 


which is impossible, since it makes F = 0. 
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3.5. The case 7 = cosh w. If H = cosh w we proceed as above to obtain, from 
(2.9) and (2.10) after slight simplification, 


tun + ty + 2tF’ = 0, (3.10) 
&+ é@r(e — 1) =0, (3.11) 
F* + (F./F), + (F,/F), = 0. (3.12) 


Let us notice by (2.15) that the surface w = 0 is a plane, so that reasoning as before 
we may assume F = ¢+ 1. 
From (2.4) we have that 


tue = 2tyt,/(t + 1) 


which may be integrated, divided by (f + 1)’ and integrated again, to give finally 


1/(¢ + 1) = a(u) + Br). (3.13) 
After division by (¢ + 1)* Eq. (3.11) becomes 
a’ + 6° +1 = Aa+t 8). (3.14) 


When differentiated with respect to u, v being constant, this leads to the linear equation 
a’(a’”’ — 1) = 0. 
We have a similar relation for 8; whence we conclude that 
a=a or (1/2)u7 + au + b, 


(3.15) 
B=A or (1/2)v? ++ Av + B. 


That both @ and 8 cannot be constant is seen by (3.13), (3.10), and (3.11); the 
equations require F = 0, which is not permissible. If @ alone is constant, moreover, the 
relation (3.14) tells us that A? + 1 = 2(@ + 8). The value of F thus obtained does not 
satisfy (3.12). Hence neither @ nor 8 is constant. After substituting (3.15) into (3.14) 
to find b + B = 1/2, making a change of variable to get a = A = 0, and using (3.13), we 
find that 


(= 2/(’+v+1)-1. 
This value of ¢, which satisfies all relations, leads to 


_ 4d? + de’) + (1+ +0*) dw’ 


dg = (2 + (cosh —ee 11 + ue bh v)} 





as the linear element when H = cosh w. Using cosh w — 1 = 2 sinh’ (w/2), writing w 
for w/2, and taking cosh? w = 1 + sinh’ w in the denominator, we find that 


dw + dv +(1+u +0) dw* 
feosh? w + (u*? + v”) sinh’ w)’ 





ds’ = (3.16) 


which assumes the proper form, we note, when w = 0. 
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The change of variable u = 7 cos 6,v = rsin @ leads to 


12 = dr’ +rdé+(i+r) dw 


(1+ (1 + 7°) sinh’ w)° 


and in this case the surfaces @ = constant are all planes, as we see by computing the second 


fundamental form. On these surfaces @ = constant we have 
. dt + dw 

ds — 7, . 2 ° , 2 
lsin’ ¢ + sinh” w 


after making the change of variable r = cot ¢. The linear element last obtained coincides 


with that for bipolar co-ordinates, whence we conclude that the spheres w = constant 
must be one of the sets of surfaces generated when bipolar co-ordinates are revolved 
about the line joining the two poles. The fact that w = 0 is plane shows that our spheres 


must be those which have their centers in the line joining the two poles. 


3.6. The case H = sin w. The final case is H = sin w, which gives 


ti+t, = F(t — 1), (3.17) 

ti. + t,, = 2F't, (3.18) 

(F,/F), + (F,/F), = F (3.19) 

when we substitute in (2.9) and (2.10), replace cos” w by 1 — sin* w, and equate to zero 
the coefficients of 1, sin w and sin* w. By (2.15) the surface w = 2/2 is a plane. Hence 
we may assume F = ¢ + 1 to obtain the result (3.13), as before. Equation (2.4) is now 


satisfied; the only ones remaining are (2.11) and those just obtained. 
. t wn) a 
By (3.13) and (3.17 
a +8,’ + 2%(at+ Bs) = 1. 


Proceeding as in the discussion of (3.15) we find that neither a’ nor @’ is zero, and that 


t= 2/(1 — uw? — v*) — 1 


so that the linear element becomes 


12 1(du°e + dv) + (1 — uw —v) dw 
as = n aa ; 
[2 + (sin w — 1)(1 ~ uw — v))° 


Continuing as for (3.16) we find that the surfaces w = constant are the other spheres 

generated by the co-ordinate lines when bipolar co-ordinates are rotated, that is, a 

set of spheres passing through a single fixed circle. This result can also be reached directly 

by inspection of ds”. 

CONCLUDING REMARKS 

4.0. Dismissal of an auxiliary condition. To complete the discussion we must show 

that (1.8) is satisfied in each of the cases considered. With this end in view we use (1.7), 
(1.5), and (1.3) to obtain, for the left member of (1.8), 

l A +A, +F A, » Aw + A, + FAQ. £1) 

{f° A’ : ere 














bo 
a 
_ 
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Replacing A by its proper value (2.3) and simplifying we find 


tu + tr. + HF’ i++ FH" 
[otetlet WY tk a ra - 
i} t+H (t+ H) 
which becomes 
tun + t.,)/8F°(t + H) — H”/4(t + A) (4.3) 
if we use (2.10). 

When H = 0, the expression in (4.3) is a function of uw and v alone, so that (1.8) is 
certainly satisfied. Similarly, if H = w” we may replace t,, + ¢,, by its value (3.2), 
whence it is seen again that (4.3) has the proper form as prescribed by (1.8). If H = 
e” the result is again true, since ¢ is zero; for H = cosh w it is a consequence of (3.10); 


for H = sin w it follows from (3.18). Thus Eq. (1.8) is satisfied automatically in all cases, 
and the co-ordinate systems hitherto obtained will actually lead to separation. What 
we have shown is that (1.8) is a consequence of the fact that the space is Euclidean. 

4.1. Non-orthogonal co-ordinates. As noted, the w-surfaces are orthogonal to the 
others, since there must be no cross derivative terms involving the variables to be 
separated. It is really not necessary, however, that the u and v surfaces be orthogonal 
to each other, although up to now we have assumed this to be the case. Thus, orthog- 
onality of the u and v surfaces is essential to our derivations, because the relations of 
Sec. 2 presuppose that the w = constant surfaces are imbedded, or at any rate can be 
imbedded, in a triply orthogonal system. Our aim now is to discard this condition of 
orthogonality as an initial assumption, and to show that the equations dependent 
thereon will be satisfied anyway as a consequence of separation. 


To this end we assume that 
ds’ = edu’ + 2f du dv + gdv” + hdw* (4.4) 


rather than (1.0) and obtain [12] in place of (1.1), 


a (gh'”  : ) (‘ h' 3 fh' ° ) 0 (8 aD 
an ( << 6.) + -1 0, — 7 0.) + aw \hi/2d 6.) = 0. (4.5) 


Upon assuming a solution RSZ as in See. 1.0 we find that 


aa. i e 8. ; 
4—-—-245,—- st + ee = 0, (4.6 

ad §S ds t d S ) 
where the terms not written involve first derivatives of the unknown functions only, 
besides e, f, g, h and R. In (4.5) and (4.6) d* is the discriminant of the u — v quadratic 
form, 


f-ag-f a7) 


Since the equation separates, there exists a function A(u, v, w) such that when 
the equation is multiplied by A, the coefficients of terms involving S are functions of u 
and v only, while the coefficients of terms with Z are functions of w only (cf. Sec. 1). 
Also the term free of unknowns must break up into a function of u, v plus a function of 
w. For our present purposes we need only the coefficients of S,, , S,. and S,, , which 
tell us that 

Ag/d’, Af/d’, Ae/d’ (4.8) 
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are functions of u and v alone. It follows that any combination of these expressions is 
also a function of u and v alone. Hence, in particular, 


g(a) _ (Af) _ a 
(49 d° \ ay) @ aii 


has this property. Multiplying (4.8) by A and using (4.9) we find that e/A, f/A and 
g/A are functions of u and v alone, so that the linear element (4.4) has the form 


ds’ = Ale,(u, v) du? + f,(u, v) du dv + g,(u, v) dv?] + hdw’. (4.10) 


Confining our attention to the terms in brackets, we see that it is always possible 
to make a change of variables, replacing u by u(u, v) and v by v(u, v), so that in the new 
variables we will have f,; = 0. This is an analytic expression of the well-known geo- 
metrical fact that every surface admits a set of parametric curves which are orthogonal. 
When such a change of variables is made, the linear element (4.10) reduces to the form 
(1.0). 

What we have shown is that there exists a change of u and v parameters which will 
make the new u and »v surfaces orthogonal to each other, if they were not originally. 
Such a change is permissible in the sense of Sec. 1.1, and may therefore be carried out 
at the beginning of the investigation. We now have a linear element of the form (1.0), 
and the foregoing derivation proceeds without further change. Thus we have completed 
the proof of the first result: 


With fixed R, if solutions R(u, v, w)S(u, v)Z(w) satisfy Laplace’s equation, and if 
separate differential equations for S and Z can be found by multiplying the equation by a 
suitably chosen function, then the co-ordinate surfaces w = constant must be orthogonal to 
the other two co-ordinate surfaces, and must consist of parallel planes, planes with a common 
line of intersection, spheres tangent at a common point, concentric spheres, the plane and 
set of spheres obtained when one set of bicircular co-ordinates is revolved about the line joining 
the poles, or a set of spheres all passing through a single fixed circle. Also, if the surfaces 
w = constant have any of these forms, and if the u and v surfaces are orthogonal to them, 
then Laplace’s equation can always be separated in the prescribed manner. 

4.2. Cases for which R = 1. Let us inquire when we may assume that &k = 1. A 
necessary and sufficient condition is that R have the form F(u, v)G(w), since in that case, 
but not otherwise, R may be absorbed in the solution S(u, v)Z(w). From the equation 
R’A = 1 we see that this condition, and hence the possibility of R = 1, is satisfied 
when, and only when, A also has the form F(u, v)G(w). We know, however, that A 
must have the form (2.3), whence we conclude that t, H, or both must be constant for 
R = 1. These considerations lead to the second result: 

We may assume that R = 1 always and only when the surfaces w = constant consist 
of parallel planes, planes with a common line of intersection, or concentric spheres. 

4.3. The wave equation. Next let us consider the possibility of separating the wave 
equation. As noted, the theory is simpler than the preceding and included in it. Writing 
the equation in the form (cf. [16]) 


Ad + ké = 0, (4.11) 


we see that F, + G, + H, of (1.1) becomes F, + G, + H, + k, and this is the only 
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change. All results previously obtained are valid here too, then, except for (1.8), which 
becomes ; 


A*(F, + G. + H, +k) = F(u, v) + H;(w). (4.12) 


But in the course of the foregoing investigation it was shown that (1.8), originally 
postulated independently as a result of the separation, is actually not an independent 
relation, but follows automatically from the others (Sec. 4.0). Though not required in 
itself for separation of the wave equation, therefore, this relation must nevertheless 
hold, in view of the other conditions. Combining it with (4.12) we see that 


A® = F,(u, v) + H,(w). (4.13) 


If ¢ of (2.3) is not constant it depends on u, say, and the same is true of F, in view 
of (1.8) and (4.13). Equating the two expressions for A, we differentiate with respect 
to u, solve for (¢ + H)*, and differentiate the result with respect to w to get 


3(¢ + H)’H’ = 0 


which shows that H is constant. Thus either H or ¢ must be constant. We therefore 
have R = 1, immediately, for the case of partial separation. In this way we obtain the 
third result: 

The wave equation separates in the sense RSZ for the three cases having R = 1, and for 
these only. 


II. Sotutions R(u, v, w)X(u) Y(v)Z(w) 


FURTHER SEPARATION OF THE EQUATION 
5.0. General. It has been supposed hitherto that the solution is only partially sepa- 
rated, being of the form R(u, v, w)S(u, v)Z(w) with a fixed function R, If we assume 
that the solution separates further to give the form R(u, v, w)X(u)Y(v)Z(w), so that 
S(u, v) = X(u)Y(v) for each function of the family, then we obtain new conditions on 
the co-ordinates. Because the cross derivative terms make separability impossible, for 
example, it is known at the outset that all three co-ordinate surfaces must now be orthog- 
onal. 
The initial stage of the separation led to an ordinary differential equation for Z(w) 
and, with m constant, to 
A’ S., A’ §S, a oa : “ 
5 a a _ ee A*F, = + A’°G, = + F,(u, v) = m (5.1) 
fs 7 s i 1g + 1g s(U, 
for S(u, v), as we see by using (1.8) and noting that each of the separated groups of 
terms must be constant. This is true because the sum must be zero, and one expression 
involves w only while the other involves u and v only. When S has the assumed form 
XY, (5.1) becomes 
a. ee 2 er Oey xi ; 
- i aw t AF, > + AG > + Fs = mM. (5.2) 
i xX + gq j + 1 Xx 1 } 4 
For separation there must be a function J(u, v) such that (5.2) separates when multi- 
plied by it. In particular the term mJ must separate, since the only other term it can 
combine with, /,J/, is independent of m. If separation is to occur for as few as two dis- 
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tinct values of m, it already implies that /,J and mJ must separate into the sum of a 
function of u only and a function of v only. In this connection it should be noted that if 
we try to pick a new J for each m, the ratio of the J values will have to be a function 
of w only and also a function of v only, in view of the condition (see below) on the co- 
efficients of X”’ and Y”’’. Thus the ratio J,/7 could depend on m alone, and the above 
considerations apply. 

5.1. The differentiated terms. It has been seen that our multiplier J(u, v) is of the 
form a(u) + B(v). When (5.2) is multiplied by this, separation must occur, and hence 
the coefficients of X”’/X and X’/X must depend on u alone, while those of Y’’/Y and 
Y’/Y depend on v alone. The condition for the second derivatives tells us that 


] i 


ds A“ }la(u) = B(v)} (dur + dv" | t dw™ } 


after a change of scale in the uw and v co-ordinates. The condition on coefficients of the 

first derivatives allows us to assume that R°A = 1, as we see by following the derivation 

of (1.7). Also we know from the first separation that A has the form (2.3), so that we 

obtain finally 

a(u) + Biv) (duu + dv’) + dw” ase 
Et, +HwP ; 


From the derivation it is clear that (5.3) is sufficient as well as necessary for separation. 

In connection with (5.3) we note that du’ and dv” have the same coefficients, just 
as in (1.5). The former result was obtained merely by introducing a change of parameters; 
it was not a consequence of separability. In the present case, on the contrary, a change 
of parameters is not permissible. The most we can do is make a change of scale, such as 
replacing u by a function of u. That the linear element has the same coefficients for 
du’ and dv” is a result which had to be proved from separability of the equation. Also 
since the u and v co-ordinates cannot be changed at will, the methods formerly used to 
simplify the equations are not available here; but there is some compensation in that 
the form of F is now known. 

5.2. The term free of unknowns. The foregoing results follow from consideration of 
the coefficients of the differentiated terms. From the other terms we obtain an equation 
analogous to (1.8), namely 


(a + B)F, = &\U) + G,(v). 5.4) 


In combination with (1.8), which we have seen is always satisfied, (5.4) gives 


a+ BA (FPF, + G. + H.) = (a + B)hy(w) + a,(u) + B,(2). 5.5) 

Using (4.3) in place of A*(F, + G, + H,) and taking F° = a + 8, we find from (5.5) 
that 

tin +t H’'(a +8, eel 

— = - a B)h, Pa! B, . (0.0) 


Si¢+H)" 4+ H) 


It may be shown that this relation, like (1.8), is a consequence of the others; we omit 
the details. The proof is closely analogous to that presented in full in See. 4.0. 

5.3. Transformations leaving equations invariant. It is convenient to note the 
changes in functions or variables which leave the equations essentially unaltered. As 
before, use of such properties permits simplification of the methods used. By inspection 
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of the linear element (5.3), or of the equations themselves, we see that the following 
transformations will not lead to any essential change, if the k; are constant: 


t — kot, (i) 
a — ka, B— kB, (ii) 
a-atk., B-—B—k., (ili) 
uf), vg), (iv) 
az2eB, Uu = ?, (v) 


DETERMINATION OF a@ AND B 
6.0. An ordinary differential equation when //’ ~ 0. Starting from (2.9), following 
through the derivation of (2.6) in detail, and using F = (a + 8)'”, we find that 


2 


a’ of id a’ a B’ i 
yy “2 fae 3% (6.1) 
2(a + 8B) 2(a + 8B) 
where c, is constant. This result assumes H’ + 0. In (6.1) we hold v constant and write 
a for a + B to obtain 











2c,a° = ala’ +c.) — (a’’ + ¢;), 


which becomes 
9 3 dp 2 
2c,a = a\p .% +c.}] — (p + ¢s) 


when we take a instead of u as the independent variable, p = a’ instead of a as unknown, 


and use a” = pdp/da. Upon introduction of a new variable S = p* = a’, this in turn 
gives 
dS 28 > 2c 
— = — ton! + dey — SS un @. (6.2) 
Qa 


da a 
The integrating factor is 1/a° and leads to 
S =a’ = 4c,a° + 2a — 3. (6.3) 
The symmetry of the equations noted in item (v) of Sec. 5.3 shows that we have a 
relation of the same type for 8. Equation (6.1) becomes a polynomial in @ and 8, as we 
= op ’ » a “ ‘ " 3 ~ 
see by computing a” and 8” from (6.3). Multiplying (6.1) by (@ + 8)” we find that 
a’ + B” must be zero when a = —8. It follows, then, that the equation corresponding 
to (6.3) for B is 
B’ = 4c,8" + 2¢.8 + 3. (6.4) 
Substitution into (6.1) leads to an identity, and hence there is no additional condition. 
Because the differential equation (6.3) also arises in the case H’ = 0, as we shall see 
below, discussion of its solutions has been deferred to Sec. 6.4. 


6.1. A partial differential equation when H/’ = 0. Turning now to the case in which 
H is constant, we absorb it into the function ¢, to obtain a simplified form of Eqs. (2.10) 
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and (2.9), Eqs. (2.4) and (2.11) remaining the same. If we take new independent variables 
@ = log ¢t and @ = log F” = log (a + 8) these relations assume a form involving deriva- 
tives only, not the unknown functions. Specifically from (2.4), (2.10), (2.11), (2.9) in 
that order we find the following: 


26.2 + 26.6, = 9.6. + O%, ; (6.5) 
du +o = but oe, (6.6) 
gun — Gre + Oi — O = Ou — 86, , (6.7) 
Bun + Ose = Wha + Po»): (6.8) 
Because of the form of @ we also have that 
ue + 0.0, = 0. (6.9) 


Let us differentiate the relation 


2bu + os) = Ou + 4,,, (6.10) 
obtained from (6.6) and (6.8), with respect to u to find that 
4(dubun + Pur) = Cae + oe . (6.11) 


Upon eliminating the second derivatives ¢,, and ¢,, by use of the other relations, one 
finds that the terms involving ¢ all cancel, leaving the equation in 6[= log (a + 8)] 


Bu. Bu + 6.5) _ Cie + Wace ° (6.12) 


It will be seen that a and 8 can be determined from this. 

6.2. Complete solution of the case H’ = @’ = 0. If 8 does not depend on v we may 
replace a by a + 8, as suggested in Sec. 5.3 item (iii), and 8 by B — B = 0. We have 
then 

O8uu = Puuu (6.13) 


from (6.12), with @ equal to log a. Letting y = 6, = a’/a, substituting and integrating 
we get 
y’°/2=wyW —c’/2 


with c constant. If c # 0 we have wy = c tan (u/2 + c’), but if c = O then 
y = —2/(u + c’). As the general solutions we thus find that 
rd »/ a 
= — [7 rae or : (6.14) 


—— - 
cos (we +c) 


In view of the equivalence (iv) of Sec. 5.3 the only essentially distinct cases are a = 
esch’u, csc’u, or 1/u’. 

6.3. An ordinary differential equation for the case H’ = 0 but a’f’ ¥ 0. We shall 
suppose now that a’s’ ~ 0, but retain the assumption that H’ = 0. With @ replaced by 
its value log (a + 8), Eq. (6.12) now gives 


a!” + B”" 


ra as = 0. (6.15) 


a’""(a + B) — 4a’’a!’ — 2a’B”’ + 3a’ 
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Proceeding as in the derivation of (6.2) we let u be constant to find that 


dS 38S Cy 
dg ~ gp t°btats=0, (6.16) 


ph? ° 7 . . . 3 . 
where s = 8’ and # has been written for a + 8. The integrating factor is 1/8° and it 
gives, after integration and multiplication by 6’, 


S =p" = A+ BB+ CA’ + Ds’. (6.17) 
By the symmetry noted in (v), Sec. 5.3, we have also that 
a” =a+t+ bat ca’ + da’, (6.18) 


where in both cases the coefficients are constant. 

Inspection of (6.15) tells us that a’ + 6’ must be zero whenever a = —8, since 
a’ and 8” are finite by (6.17) and (6.18). If a” is given by (6.18), then 6’" is uniquely 
determined as 


B" = —a+ bB — cB’ + ds’. (6.19) 


It may be verified that (6.15) is now satisfied, and hence that there is no additional 
condition. Since these equations are equivalent to (6.3) and (6.4), obtained for the case 
H' = 0 (cf. Sec. 6.4), their solution completes the determination of a and 8. It is seen, 
incidentally, that if a(u) satisfies (6.18), then the function —a(iv) will satisfy (6.19). 

6.4. Canonical forms of the equation. Let us make use of Sec. 5.3 to simplify (6.18) 
and (6.19). Writing —a@ for a will make d > 0, if originally d < 0; and replacing a by 
Aa + B, u by u/C, we get 


A’®a’'C? = a+ b(Aa + B) + c(Aa + B)*? + d(Aa + B)*. (6.20) 

Supposing that d has been made greater than zero, we set 
a+ bB + cB’ + dB* = 0, (6.21) 
b + 2cB + 3dB* = dA’, (6.22) 


and choose C = (Ad)'” to obtain the canonical form 
a” =a +o’ +a. (6.23) 


The corresponding substitution for 8 must be 8 — AS — B to be permissible; it gives 
(6.23) with a replaced by 6 and A by —A. . 

We see that A, B, and C may be assumed real (this is important), and that A = 0 
is necessary only when there js a triple root of the original equation. These observations 
follow from the facts that the left-hand side of (6.22) is the derivative of (6.21) and that 
d > 0. For a triple root the canonical form is 


a’ = a’, B” = B*. (6.24) 
If we suppose now that d = 0 but c + 0 we find, in a similar way, that 
a =a+dre+1, 6 = —-8+284+1. (6.25) 


This form includes the double root case, which can be shown, however, to be impossible 
anyway. One would at first expect a + sign in front of a@’’, but this is accounted for by 
interchanging a and £. 
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Next c = d = 0 gives 
a =a, pf =8B (6.26) 
and a’ = ~’ = Ois the last possibility. In Eq. (6.23) we must distinguish the three cases 
|X| < 2,{A| = 2,|A| > 2, whereas in (6.25) we have only the one case || > 2, 
since a’ and #’ must be real simultaneously. 

6.5. Solution of the ordinary differential equations for a and 8. We have seen that 
whenever H’ or a’8’ + 0, the functions a and 8 must be solutions of one or another of 
the canonical differential equations in Sec. 6.4. These equations, though non-linear, may 
all be solved by elementary methods. The work is rather tedious, particularly since one 
must be careful to keep all possible solutions. We, therefore, content ourselves with a 
single example. From (6.23) we find that 


du = fala” + rA¥a + 1)] 7? da 


which leads (among other possible expressions) to [14] 


/ 4\1/2 Ty ae ay 
utc= e t) f E (7 — r’)”’”.” 6 | sin’ ¢@ = o—, 


when the roots 7, 7’, 0 are real and unequal, that is, when | A | > 2. Here F(z, p) is the 
elliptic function of the first kind. Since we want a@ as a function of u rather than the 
converse, we introduce the Jacobi elliptic functions to find that 

” ol > 11/2 

7 = sn" (ar' “(u +c), aT = Fr ) 

ae 
with a similar expression for 8. For r and 7’ we substitute the proper values in terms 
of \ to obtain, with a new wu, a one-parameter family. This is for the case0 < 7’ < a < 1. 
Other cases are similarly. treated, and lead to two additional expressions when | \ | > 2. 
The cases | A | < 2, which gives conjugate complex roots, and \ = 2, which gives equal 
roots, are dealt with in the same way; the latter and all others in (6.24)—(6.26) give 


elementary functions 
DETERMINATION OF f 


7.0. The partial differential equations. Up to this point we have shown that the 
linear element of orthogonal co-ordinates in Euclidean space must be of the form (5.3) 
if Laplace’s equation is to separate. In addition H is one of the values (2.7) and @ and 
B are each one of the values obtained in Secs. 6.0—6.5. It remains only to determine 
tin (5.3), and we proceed to this question forthwith. 

Treating first the case H’ = 0, our aim is to compute the expression involving F 
on the right of (2.9). Assuming that a’B’ or H’ ¥ 0 we have (6.18) and (6.19). Differ- 


entiating both sides with respect to u and dividing by 2a’ we get 


9) 
” 0 2 c 4 
a = aa + ba + ms (7.1) 
: ae , 1/2 ; 
with a similar expression for 6’. Since F = (a + £8)” we have, for the expression de- 


sired, 


F \ 'F \ i l a’! + g** a! 2 Br 
F a, (F } 7 Al atB  (@ + 5 
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which reduces to 


) 4(B) «a+o8 a9 
in view of (6.18), (6.19) and (7.1). Here a is the coefficient of a* in (6.18), and hence 

| for the cases with canonical equations (6.23), (6.24) but a = 0 for (6.25), (6.26). 
The above assumes (6.18), which is valid only when H’ or a’s’ # 0. Assuming now 
that H’ = ,’ 0 and writing a for a + 8 we find that the expression on the left of 
7.2) becomes (a’/a),/2. For the case H’ = B’ = 0 it has been shown that a must be 
one of the three functions mentioned at the end of Sec. 6.2. By direct calculation we 
a every time, and hence (7.2) is valid with a = 4. Thus the ex- 


a 


find that (a’/a),,/2 
pression is known in all cases. 


When H = 0 we find that 


t f la , 
te + thy = A + 7 (a + £) (7.3) 
l — 
ott, +2-— (7.4) 
from (2.9) and (2.10), af making use of (7.2) and replacing F* by a + 8. These re- 
lations tell us that 4; + ¢ = 0 whenever a 0, so that in every non-cubie ease if H is 
constant ¢ must be constant also. The only remaining relations are (2.4) and (2.11), 
which become 
21 _ &itat 7.5 
— at+8B / (75) 
a’l, — B't, 
t.. - t,, = a 7. 
- a+86 (@ 6) 
for the present situation F* = a + 8B. Equations (7.5) and (7.6) are valid for all H. 
When H is not zero the corresponding form of (2.9) and of (2.10) is obtained by 
simply setting F a + 8 in the relations of Secs. 3.1—3.6. Specifically, if H = w* we 
have that 
ti + t., = 4(a + 8B), (7.7) 
(t.. +t.) E+E; (7.8) 
if H e” then ¢ 0; if H = cosh w, then 
t., + t., + 2ta + 6) = 0, (7.9) 
éj@+@+(at+ pA — 1) =0; (7.10) 
and if H = sin w, then 
ti+t=(a+ ay — I), (7.11) 
uu + t., = 2(a + Be. (7.12) 
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It may be verified that both (2.9) and (2.10) are satisfied identically when, for a 
given H, the function ¢ satisfies the appropriate pair of equations from (7.3), (7.4) or 
(7.7)—(7.12). These relations plus (7.5) and (7.6), then, give the necessary and sufficient 
condition for separation in Euclidean space. 

7.1. An indirect method of solving the equations for ¢. ‘lo obtain all possible linear 
elements we must have all possible values of ¢; in other words we must have the general 
solution of the above differential equations. This general solution is not easy to obtain 
directly and hence we use an indirect procedure, which depends on the fact that a and 
B are already known for every case. 

Given a particular linear element ds“ we note from (5.3) that all others using the 
same a and £ are obtained from this one by means of a conformal transformation, with 
ds’ now regarded as the metric of Riemannian V; . Thus, since a and 6 are to be the 
same in both cases we can change ¢ and H only. The two elements are, consequently, 
proportional. By a theorem of Liouville [15], the most general conformal transformation 
which preserves the Euclidean character of the space—i.e., preserves the relations of 
Sec. 2—is a rigid motion, a reflection, or an inversion. Only the latter is of interest here. 
Hence given a and 8, if we can somehow discover just one admissible ¢ and H, then all 
others can be found from this by inversion.* 

A single value of ¢ is found without much trouble from the equations; in many cases 
it is obtained by inspection. Or one may refer to [3] and [11], which give examples of 
Euclidean linear elements with our a and 8 for each case. Knowing the single solution, 
we get the others by inversion as outlined above, and it remains only to see whether all 
inversions are permissible. : 


CONCLUDING REMARKS 


8.0. Concerning inversion. It was found by actual trial that all inversions of per- 
missible co-ordinates were themselves permissible, if by permissible we mean that the 
space is Euclidean and the equation separates. The new linear element had the form 
(5.3) and the other conditions (7.3)-(7.12) were always satisfied, independently of the 
center of inversion. We propose now to give a direct proof that this must necessarily be 
the case. 

First, since the original linear element was permissible, it had the form (1.5) with 


F? = a+ 8B. Hence, because an inversion is a conformal transformation, the new linear 
element will also have this form with a different A, say A. The equation will separate 
if R, the new value of R, is given by R°A = 1. This is true provided only that the rela- 


tions (1.8) and (5.5) continue to hold in the new system, and we shall see that such is 
the case. 

Since the original linear element was permissible, it was obtained from a co-ordinate 
system in Euclidean space. Consequently the new one, having been found by inversion, 
will have the same property. Now the relations (1.8) and (5.5), which are the only 
remaining conditions for separation of the equation, have been proved to be satisfied 
automatically whenever the space is Euclidean. They are therefore satisfied by the 
new linear element, and hence the equation separates. Moreover the fact that A has 
the form (2.3) was also deduced from the Euclidean relations, and hence will persist 
for A; and the same is true of all the differential equations for ¢ and H. Thus the linear 


4A similar use of Liouville’s theorem is made in [8]. 
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element obtained by inversion satisfies every one of our conditions if the original linear 
element does. These remarks, which apply to partial as well as complete separation, 
complete the proof of the fourth result: 

If the above procedure can be applied to the resulting differential equation for S to get 
further separation, so that the solutions have the form R(u, v, w)X(u)Y(v)Z(w), then the 
co-ordinates must all be orthogonal and must be either toroidal co-ordinates, or the well- 
known cases (I—III but not IV in [3]) giving separation of this type with R = 1, or else the 
co-ordinates obtained by inversion of these in a sphere. 

We have also proved that: 

The set of co-ordinates giving separation RSZ or RX YZ, as the case may be, is closed 
under inversion. 

In other words, if a given co-ordinate system is in the set, then the system obtained 
by inverting it in any sphere will likewise be in the set. This behavior of course is not 
found when we confine ourselves to the case R = 1. 

It must be mentioned that these results on inversion, though interesting mathe- 
matically, are of slight practical importance. The way one would actually solve Laplace’s 
equation in co-ordinates which are inversions of standard cases would be to solve the 
standard case and then invert the solution. The fact that the equation could be solved 
directly by separation of variables in the new co-ordinate system, though true, would 
not be used in practice. 

8.1. Conditions for which # = 1. We shall have R = 1 when, and only when, it can 
be absorbed in the product XYZ, that is, when R = p(u)q(v)r(w). In view of the condi- 
tion R°A = 1 in Eq. (1.7) we see that A too must have this form. Such a condition, 
when combined with (2.3), makes H constant and F of the form r(u)s(v), or else it makes 
t constant. It is found that these conditions hold only in cases which occur in [3]. Thus, 

Laplace’s equation separates in the sense XYZ when, and only when, the wave equation 
so separates. 

The ‘‘when”’ part of this result is well-known, but the ‘only when’’ part is believed 
to be new. 

8.2. The wave equation. In Sec. 4.3 we considered the conditions under which the 
wave equation would separate partially, to give solutions of the type RSZ. Turning 
now to the case in which there is complete separation RX YZ, we encounter a difficulty 
when we try to show that F* = a(u) + B(v). The previous argument depended on the 
fact that one term of (5.2) involved the separation constant m while the other did not. 
For the wave equation, however, we may have k depending on m, and this method 
therefore cannot be used. Confining our attention to the case H = 0, which we have 
seen is the only case in which F need not be constant in the wave equation, we shall 
re-examine the separation procedure.° 


Since Sec. 5.1 makes J = F? = G’, the same equations are found, in general, with F? 
replacing a + 8. But (5.4) becomes 
F°\ A°(F, + G. + H. +k) — m} = a2(u) + 8.(v) (8.1) 


when F, is replaced by the F’; of (4.12). Observe that H = 0 in (2.3) makes (1.3) depend 


’This discussion was revised in proof (April, 1949); the earlier version was incorrect unless we had 
kim, = kom, in (8.2). 
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on uv and v alone, by (1.5) and (1.7). Hence H; of (4.12) is constant, and may be absorbed 


in F;. 
If there is separation for three pairs of values (m; , k;) such that 


] 


l m hi; 
1 m x () 8.2) 
l m / 


then the three equations (8.1) may be solved for F”, whence we obtain the desired result 
| : 


F° =a + £. 
When the determinant (8.2 
complete result has not been obtained. It appears that separation can occur for co-ordi- 


is zero, which can happen in practice (ef. Ref. [16]), the 


nates not obviously reducible for the known cases, and that J = a + @ is no longer essen- 
tial. Also a rather long investigation shows that 6 = log F” must satisfy 0,, + 0,, = cF’ 
und that the constant c is zero if any minor like m,k. — mk, vanishes. 


“TI” Bes 


when F° ¥ a 





This latter condition gives ¢ 
which the determinant 


Q by (6.6) and (6.8 


, whence we conclude that F” = a + 8B. 


(8.2) vanishes, but no minor does, 


IS 





The general problem in 
reserved for later investigation. 

Now that we have obtained the proper form for F we can use all the results of the 
preceding sections. These results, which certainly represent a necessary condition on the 
co-ordinates for separation of the wave equation, tell us that we have ¢ or H constant 
only in the well-known cases giving R = 1. Consequently, 

The wave equation separates in the sense RX YZ only in the known cases for which it 


separates in the sense XYZ, provided (8.2) holds. 
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LINEAR EQUATION SOLVERS* 


BY 
F. J. MURRAY 


Columbia University 


1. There are two types of devices for the solution of simultaneous linear equations 
which have been developed. Suppose the given system of equations is 


>> a,;2; = b;, i=1,2,---,n. (1) 
?=1 


In one type the b; are fed into the machine in such a way as to drive the unknowns 
zx; to their correct values. In the second type, the x; are not driven by power supplied 
from the constant inputs but reach an equilibrium situation corresponding to the 
solution by a process of adjustment. We are concerned in this article with the operating 
conditions for this last type of machine when the adjusting process is determined by 
a linear operator with constant coefficients. 

We suppose that each a,; can run independently through a real range which contains 
the origin. Thus, the determinant may be zero and any matrix can be represented by 
a suitable choice of scale for the unknowns and the b,; . Adjusting machines which are 
stable even when the determinant is zero may be designed. For instance, a block dia- 
gram is given in the author’s book’ for such a machine. Another example is the set-up 
described by Goldberg and Brown’ which will insure stability when a certain type of 
feedback is used. 

However, in each case the coefficient network is duplicated. In the present article, 
we point out that if an adjusting type of machine is to operate successfully whenever 
the determinant A is not zero, then the square of the determinant must enter the indicial 
equation of the equations of motion for the machine. This necessary condition for 
successful operation rules out any linear feedback which does not involve using the 
a,, twice. This result generalizes certain aspects of the necessity argument indicated in 
Goldberg and Brown. 

In Sees. 2 and 3 below, we describe precisely the type of machine we are concerned 
with. These machines may function continuously or in discrete steps. In Sees. 4 and 
5 we obtain necessary and sufficient conditions that the machine should operate suc- 
cessfully in all cases where a solution is uniquely determined. These conditions are 
analogous to stability conditions for a linear network. In the case of a continuous ma- 
chine, this analogy is readily established; in the case of a discrete step machine, these 
ope rational conditions are obtained by considering certain parts of the theory of linear 
difference equations. In See. 6 we prove the mathematical theorem upon which our 
result is based. It is shown in See. 7 that an adjusting machine with a linear feedback 
network, which is independent of the coefficients of the equations, will not always 
operate successfully. Section 8 contains the mathematical basis for Sec. 5 which is 
concerned with discrete step machines, 

*Received May 28, 1948. 

F. J. Murray, Mathematical machines, Kings Crown Press, New York, 1947 p. 92 (1st ed.), pp. ITI 


20-21 (2nd ed 
Ii, A. Goldberg and G. W. Brown, J. Appl. Phys. 19, 339-344 (1948). 
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2. A mathematical machine can be regarded as a combination of computing com- 
ponents, each of which performs a specific mathematical operation. Each component 
has various inputs and a specific output. A combination of components can be used to 
evaluate a formula or function of the inputs. In a machine for solving a system of linear 


equations 


2, a;,t%; + 6; = 0 (1) 


the coefficients a;; and the constants b, are inputs whose values do not vary during 
the operation of the machine. 

In an adjusting type machine for solving linear equations, there are variables to 
represent the xz; . Each of these is associated with a unit of the following character. 
Each unit has an input X; and an output x; . The output z; is suitable for use as an 
input in the computing components of the machine and the input X,; must correspond 
to an output of computing components. The relation between the input X, and the 
output z, is an operational one, L,(z,) = X;, . Various possibilities for L; are discussed 
below; L, may depend on 7, it will always be linear. 

The adjusting type device functions as follows. While the xz; units are inactivated, 
the a;, and the b; are entered into the machine. Presumably at this point the values 
of the xz; do not constitute a solution of the given system of equations (1). Now, however, 
the x; units are activated. Various combinations of components compute the errors in 


each equation 


ée- = > aya; + dD; (2) 
7=1 
and these, in turn, are used to compute the X, as some function f(e, , --- , €,) of the 


errors. The inputs X, of the z, units cause the latter to vary so as to approach the 
solution of the system (1). Since as we have pointed out in Sec. 1, for this purpose 
certain f values exist which are linear combinations of the e; , we will assume that f is 
linear in the ¢; : 


X;= Dd kise; (3) 


where the k,; are constants in the sense that they do not change during the process in 
which the zx, are adjusted to the correct value. 
3. We have not specified the nature of the operators L, , 


L(x) = X;, (4) 


or the method of functioning for the components. In general, there are two ways or 
manners in which a machine of this type may operate: (A) The adjusting process pro- 
ceeds continuously. Each component has continuous inputs and output and the L, 
are differential operators. We suppose the coefficients in each L, are constants and 
that the coefficient of the highest derivative is 1. (B) The adjusting process proceeds 
in discrete steps. L, is a linear difference operator with constant coefficients but the 
components may have discrete inputs or continuous inputs. Again the coefficient of 
the highest order difference in J, is 1. 

The two examples cited in the introduction operate in the (A) manner. An electronic 


digital computer programmed to operate in the sequence indicated by Eqs. (2), (3), 
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(4) would operate like (B). A non-essential generalization permits us to include the 
manually adjusted type of equation solver such as those described in the author’s book 
Mathematical machines (2nd edition, pp. III 16-20; Ist ed., pp. 87-91). A full cycle of 
the adjusting process in these machines corresponds to one step in the sense of (B) 
above. However, the adjusting equations (4) are to be replaced by a set of relations in 
the form 


L(x) = X,, 
L.(a, ’ X2) = X2 ’ (4’) 


L(x; doe » &D _ 7. 


where each L is linear in the x; . It will be seen that the generalization represented by 
(4’) does not affect our argument. One may mention that Z, may be taken simply as 
the differencing operation A. In general L; contains x; only in the form Az; so that 
the first equation determines Az, , the second Az, and so forth. 

We are interested then in the type of machine represented by the sequence of 
equations (2), (3), (4) or (4’). The obvious problem that appears here is concerned 
with the choice of the k,; in (3) and the LZ, in (4) so that the machine will work, i.e. so 
that the machine will adjust itself to a solution. As we have pointed out, sufficient 
conditions for adequate operation in all circumstances are known, but these require 
that the k;; depend on the a;; . This means that the a,; must be used in computing 
e, and also again in forming the X; . 

Our objective is to establish that this double use of the a,; is necessary if the machine 
is to function in all cases in which there is a solution. We suppose that the a,; are per- 
mitted to assume independently all values in an interval which includes zero, say, for 
instance, from minus 1 to plus 1. If this is true then, by suitable choice of scale for the 
equations and the unknowns, any system of equations can be represented in the machine. 

The above statement seems to neglect the case in which the a;; appear digitally or 
as decimal fractions. However, when we are given a device in which the a;; vary dis- 
cretely in the components, we can regard each component as replaced by a continuous 
component of perfect accuracy and apply our argument to it. Now suppose in the 
idealized machine, we find a region of non-operation or instability. Then, in general, 
the original discrete machine will have permissible values for the a,;; which fall in this 
region and, of course, it will also be unstable. 

Let us end this section by pointing out certain conclusions concerning the above 
mathematical setup which one may reach from the assumption that the machine will 
operate successfully as an adjusting device. 

1) If the machine operates in the manner (A) and the ZL, are differential operators 
such that 


Lz) = 2{° + Laie? + e+ + La, 


then /,; = 0. 
2) If the machine operates in manner (B) and 


L,(x;) = A‘z; + 1, ,A°~ "2; + om + LX; ’ 
then l,; = 0. 
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From Eqs. (3) and (4), in either case, we get 
Lila) = >> kije; - (5) 


Now suppose we are solving a system of equations in which 2; is not zero. Then after 
an adequate time interval in a successfully operating machine all the e,; will be small, 
zx; will be close to its true value, all the differencing operators A‘ and differentiating 
operators (d/dt)" will yield a small result, and these equations will be inconsistent with 
the assumption that 1,; # 0. The case (4’) is readily treated on an inductive basis. 

4. Consider the Eqs. (5) in the case (A): 


1 


L(x) = 2° Hee H+ a ti 


= 2, hue (5A) 


wad Ban. 44) 


Equation (5A) can be written in the form 


L(z;) — > ( z i308; Je = p k,;b; = (say) B, 


The operators L,; can be treated as numbers, and we may employ the process by which 
Cramer’s Rule is usually established to eliminate all but one unknown z. Consequently 
each z satisfies a differential equation in the form 

Vz, = C, (6) 
where V is a differential operator, 


Vi, =2;" + Diz ob eee + Diz, , 


and each C; is a constant which of course depends on the b; . The coefficients D do not 
depend on 7. Clearly D,, is the determinant of the matrix with elements b kiiQja- 
(This follows from the fact established in the previous section that l,; is zero.) Thus 


D,, = KA, (7) 


where K and A are, respectively, the determinants of the matrices with elements k;,; 
and a,; . 

We proceed next to obtain the condition mentioned in Sec. 1 for successful operation 
in case (A). This condition is concerned with the algebraic equation 


wu” + Dw" + --- + D, = 0, (8) 
which is usually referred to as the indicial equation of the homogeneous differential 
equation Vz; = 0. 

We first obtain the general solution of Vz, = 0. Suppose yu, , --- , w., are the real 
roots of (8) and that yu; has multiplicity 7; . Also suppose that a, + 7B, , a, — 78;, °°: , 


a,, + 78,, a,, — 78,, are the complex roots of (8) and suppose that a; + 78; has multi- 
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plicity u; . Then it is well-known that the general solution 2¥ of Vz, = 0 can be written 


in the form 
zt = > 2 M,.,t’~* exp (ut) + 2. + t’~' exp (a,t)(P;,, cos B.t + Q;,, sin Bt). (QA) 


On the other hand if s» is such that D,,_,, is the D with highest subscript which does 
not vanish, then a particular solution of (6) is obtainable in the form X,,t’’ and the 
general solution of (6) can be written in the form 


t’? + o*. (10) 


Now suppose that the system of equations has a unique solution in which no unknown 
Y,; is zero. For the machine to function correctly under these circumstances, each 2, 
must approach the correct constant value Y,; as ¢ approaches infinity. Now one can 
show that an expression of the type (10) will approach a constant value as ¢ approaches 
infinity if, and only if, s) = O and all the «4, and a, are negative, i.e. if (8) is stable. 

Thus we have shown that: 

An adjusting machine of type (A) will operate successfully if, and only if, all the roots 
of the indicial equation (8) lie in that half of the complex plane for which the real part of a 
number is negative, whenever the system of equations ts non-singular. 

5. We now wish to go through the analogous discussion for case (B). In this case, 

the L 


are difference operators 


L(z;) = A‘r,; + l, -_ "Xt; + itil + Ras ;Ax; ’ 


and we obtain 


n 


LAr) = >> kise; . (5B) 


As before, we also obtain by elimination 
Vz; =C,, (6B) 
where 
Va, = A”r, + DA" x, + +++ + Daa, . 


The statements concerning D,; and D,, are the same as before and, in particular, (7) 
holds. 

The customary method of handling a homogeneous difference equation Vz; = 0 
with constant coefficients is to look for solutions in the form 2z,(t) = a‘. For these we 
have Az, = a‘*' — a’ = (a — I)a‘ = (a — 1)z,(t) = ua,(t) for up = a — 1. Thus if 
we permit « to have this meaning, we can obtain the indicial equation 


w+ Dw™' +--+ + D, = 0. — 


If we had two machines, one of which operated in manner (A) and the other in manner 
(B) with, however, the L,; and X, being such that the coefficients L,;; and k;; were the 
same in each case, then the Eq. (8) would also be the same for each case. 
Corresponding to every distinct solution of (8B) we have a solution a‘, where a = 
u + 1. For multiple roots the same technique as that of ordinary differential equations 
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permits one to find r linearly independent solutions if the root has multiplicity r. For 
instance, if uo has multiplicity 3 then there is a polynomial ¢ such that 


we” + Dy" + +++ + Dy = (u — wo)*o(u) 
and 
Vu + 1) = (uw — wo)*o(u)(u + 1). 


From this we can infer that 


V(uo + 1)‘ = 0 


2 196+ DL... = 0, 
Ou 


a” in 
ae [V(u + 1)'],-.. = 0. 
Since A and 0/du are commutative operations, we can infer that 


v(2 l(u + 1'Teens) = 0 
Ou 


and 
u(% [(u + 1)'Jeen = 0. 
7m 


Thus (uo + 1)‘, (uo + 1)‘ " and t(t — 1)(uo + 1)‘ are solutions of Vz; = 0 when 
Ho is a triple root of the indicial equation. Since linear combinations of solutions are 
also solutions (uo + 1)‘, t(uo + 1)‘ and #’(u. + 1)‘ are equivalent to the given set if 
Ho + 1 does not equal zero. If uw + 1 = O, then xz; = 1 is a solution; if uw = —lisa 
root of multiplicity 3 of (8), then, of course, the analogous equation for a has zero as 
a triple root and thus x; = 1, ¢ and @ are solutions. If yu» is a complex root and the D, 
are real then ji is also a root and solutions are obtained by taking the real and imaginary 
parts of the above solutions. 

It is, of course, characteristic of the theory of linear difference equations that it 
parallels the theory of linear differential equations except that constants are replaced 
by periodic functions in the solutions. (There is another departure indicated below.) 
However, we need only consider the variable ¢ to have discrete values 0, 1, 2, --- ; 
then periodic functions need not enter our discussion. For this situation we give a 
sequence of lemmas which yield results analogous to those we have used in the differ- 
ential equation case. These are proved in Sec. 8 below. 


Lemma A. Let 
Vz=0 (6B’) 


be a linear difference equation with constant coefficients. Let the variable t take on only the 
values 0, 1, 2, --- so that a solution of (6B’) is an infinite sequence 


{x(0), z(1), x(2), Pare }. 


There are m linearly independent real solutions of (6B’). 
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Coro.uary. For the solutions obtained in Lemma A, we can suppose that the first de- 
terminant ts not zero. 


Lemma B. Under the hypothesis of Lemma A, every solution of (6B’) is a linear combi- 
nation of the m linearly independent solutions obtained in Lemma A. 


Lemma B specifies the solution of the homogeneous equation Vz; = 0; we can 
readily show that this is in a form analogous to (9A). For to each solution yu of (8) we 
can find a principal value of the logarithm such that 


1 + uw = exp (a’ + 7p’) 
or if 1 + uw is real and positive, we can write simply 
1 + » = exp (y’). 


If s) is again such that D,,_,, is the D with highest subscript which is not zero, then a 
particular solution of (6B) is readily obtained in the form 


t! 


=e @ 1 oon. 
, *** Sol(t — 8)! 


The argument used in the previous section now shows that if the machine is to 
function satisfactorily, then when the solution is unique one must have all y’ and a’ 
negative. However, we are now one step removed from Eq. (8) and the condition for 
satisfactory operation of the machine is now that all the roots of (8) should be in the 
unit circle with center at 1 = —1. Thus we have proved: 

An adjusting machine which operates in the manner (B) will operate successfully tf, 
and only if, all the roots u of the indicial equation (8) lie in the unit circle with center at 
uw = —1, whenever the determinant of a system (A) is not zero. 

6. We have established two facts which hold for adjusting machines quite irre- 
spective of whether they operate in manner (A) or (B). 

a) The last coefficient D,, of the indicial equation is in the form KA where A is the 
determinant of the given system of equations. 

b) If any root u of Eq. (8) has a positive real part when the given system of equations 
is non-singular, then the machine will not operate successfully for this system. 

We now wish to establish the following theorem upon which our conclusions are 


based. 


Turorem. Let Eq. (8) be the indicial equation of the Eq. (6) obtained by eliminating 
all but one x; from the Eqs. (4). Furthermore, let us suppose that tf the system of equations 
(1) is non-singular, then none of the roots » of the indicial equation (8) has positive real 
parts. Then D is divisible by A*, where A is the determinant of the system of equations (1). 


Proor. Let us suppose that our system of equations has been chosen in the following 
way. Let A(A) be the characteristic determinant of the matrix and suppose we choose 
a set of a,; so that all the latent roots, i.e. the roots of A(A) = 0, are real and distinct. 
A symmetric real matrix of this type is readily found. Furthermore, for real a;; one can 
show that for small changes in the a,; the roots remain real and distinct. This is readily 
seen from the graph of A(A) which will have n — 1 maxima and minima whose ordinates 
alternate in sign if, and only if, A(A) = 0 has n distinct roots. Thus we have an open 
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region in the n° dimensional a;; space in which the characteristic roots are real and 
distinct. 
Let us, however, for the moment regard all the a;; as fixed but consider all systems 


of equations in the form 


a, — Aja, + --- + a,,2, + b, = 0, 
(1.A) 
Gnjt,) +--+ + (a, — Ax, + OD, = O. 
Equation (8) can now be written in the form 
wu” + D,(A)u" | + ee + *<D, (A) = O, 
where now the D,(A) are polynomials in a,; and \. We know that D,,(A) = KA(A). 


(K is constant, but only in the sense that it does not change during the operation of 
the machine, i.e. in the sense a;; and X are constants.) 

Now suppose }’ is a root of A(A) = 0. Then since D,,(A) = KA(A), uw = 0 is a root 
of the indicial equation for the system of equations (1.A). Let p be the last integer such 
that D,,_,(d’) does not equal zero. Since D,,-, is a polynomial in \, Dz -.(A) = 
(A — X’)D,,-.(A) for k < p where D,,_,(A) is also a polynomial in X. 

Now suppose D,,(\’) # 0. Let us make the usual construction of the Newton polygon 
in order to obtain a series expansion® of the roots yu of the indicial equation in powers 
fi — WV: 


w= e,(A — N’)° + ex(A — 0) + 


We readily find that (0, 1) and (p, 0) are vertices of the polygon. Let 6, , --- , @, denote 
the p roots of the equation 


DAY). - 6D. AN) = 0. 
We have for each k = 1, --- , pa solution in the form 
up = 0,(X — r’)'? + 


Now if p = 1, 4, is real and 6,(A — 2’) changes sign as A varies through \’. This would 
yield a uw with a positive real part when A(A) is not zero. But this violates our hypoth- 
esis. For p = 2, we can choose \ so that both 6,(A — d’)'” and @,(A — X’)'” are real 
and one will be positive. This will again destroy the stability. For p = 3, 4, --- we 
can take \ — 2’ positive and since at least one 6, has a positive real part, the real part 
of uw will be positive, and again we have contradicted the hypothesis. 

Thus the assumption that D,,(\’) # 0 contradicts the hypothesis that the roots 
have negative real parts when A is not zero. Hence D,,(A’) = 0 and D,, has a factor 
(A — d’)’*. 


Now let us divide D,,(A) by A*(A), considering them as polynomials in \ and the 
a,; . Since the coefficient of \*” in A*(A) is one, we find that there are two polynomials 


@ and R such that 


3Cf., for instance, K. W. S. Hensel and G. Landsberg, Theorie der Algebraischen Funktionen Einer 


Variabeln, Teubner, Leipzig, 1902, pp. 39-52. 
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D,() — Q(A)A?(A) = RA), 


where R(A) has degree less than 2n in X. 

This holds for all matrices {a,;} in the {a,;} region we decided upon, within which 
all the roots \’ are real and distinct. However, in this region for each choice of {a;;}, 
R(A) is divisible by (A — 2’)* for every root \’ of A(A) = 0, and since the degree of 
R(A) is less than 2n, R(A) must be identically zero in for {a,;} fixed. Since this holds 
for every matrix {a;;} of the region, R must be zero as a polynomial in a;; and \. This 
establishes the theorem. 

7. The above necessary criterion for successful operation permits us to show that 
a simple adjusting machine, with a linear feedback in which the a;; are not used, cannot 
function successfully for all systems for which the determinant is not zero. For here 
D,, = KA where K does not depend on the a;; and thus D,, is not divisible by A’. 


8. Lemma A. Let 
Vi =0 (6B’) 
be a linear difference equation with constant coefficients. Let the variable t take on only the 
values 0, 1, 2, --- so that a solution of (6B’) is an infinite sequence 
{x(O), x(1), 2(2), ---}. 
There are m linearly independent real solutions of (6B’). 


Let us first ignore the restriction that the solutions be real. The above argument 
indicates how m solutions can be obtained by means of the indicial equation. We must 
establish then the linear independence of the m solutions obtained: 


1 


f2"O), = 2° (1), +}, 


fr'”(0), x(1), «°°, 


is oe, raed (| a 


tegarding these solutions as an infinite matrix, we must show that the rank is m. If 
we do this, a finite submatrix must have rank m and have linearly independent rows. 
Consequently, the full matrix has linearly independent rows. 

If the roots of (8) are all distinct then, of course, the determinant of the first m 
columns is the well-known Vandermonde determinant 


uw +1 vee (wu, + 1)"" 
l Mo + 1 tae (uo + 1)""" 

= [] (us — wy), 
1 tm + I ee (Hy 





which is not zero. 
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The case of multiple roots is' treated by induction. We suppose that yo is a triple 


root and that the first three solutions correspond to yo . In addition, for the moment 


we suppose that all other roots are distinct. Then the first three rows of the matrix are 


| Ho + 1 (uo + 1)° (uo + 1)° (uo + i)* 


40 1 2(uo + 1) 3(uo + 1)’ 4(u + 1)° tee 
lo 0 2-1 tj. +D 43649 --« J, 
In the distinct root case, we have 
1 0 wotl tes (wy +1)" | 
l Hi + 1 sets (ui + | fin | 
] Me + 1 se (u2 + 1)" | 
| 


= (uo — Hi)(o — oe) — we) TT (wo — ws) TI] (ns — ws) TT Ge — wd) T]  —- 2) 
where 7, 7 do not equal 0, 1, 2. If we let u; = uo + h, the elements in the second row of 
the determinant become 
(ur + 1)’ = (wo + 1)’ + AlG(uo + 1)" + 1/2755 — (wo +1) 7] +--+» + hi. (11) 
Now we can subtract the first row of the determinant from this second row, divide both 
sides by h and take the limit as h approaches zero. Then we obtain 


| ] Ko + ] (uo +> 1)’ eee (uo + re | 

| 
| 0 l 2( uo > 1) coe (m a, 1) (uo “4. ta 
| | 
| ; | 
| 1 fe + I (uz + 1)° nei (us + ai | 


| 


| 
' 


= (uo — we)” [] (uo — ws)? T] (ue — ws)? I] (us — a). 


We next let uw. = wo + A and the third row is in the form (11). Then we subtract 
the first row and fA times the second row from this third row, divide both sides by h’, 
and take the limit as h approaches zero. Thus we get 


11) wotl (wo+1) -:: (uo + 1)" | 
‘Oo 1 eS , ae (m — 1)(uo + 1)""? | 
1 } 
“10 #0 2-1 s+ (m — 1)(m = 2)(u + 1)""* | 


= I] (ito = hi) I] (us = Mj). 
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Since the determinant has been expressed as a product of differences of distinct numbers, 
it is not zero. It should be clear how the cases of higher multiplicity are treated and also 
how one could successively deal with the case of two multiple roots, three multiple 


roots and so forth. 
Thus, in every case we have that the solutions previously obtained are linearly 
independent. If we have complex roots, we write the determinant in the form 


}1 owttl ov G+ 


1 opgtilooes. = @+1)"" 





| 


i.e. every line containing a complex number is followed by the corresponding line for 
the complex conjugate. Now a determinant in the form 
A, Az As A, 


i is fk 


B, B, B; B, 








B, B, B; B, 
can be converted into one whose first row consists of the real parts of the A, and the 


second row the imaginary parts, and whose third and fourth rows depend in the corre- 
sponding manner on the B, , by multiplying on the left by the determinant 








1 1 
5 9 0 0 
—1 t 
= 9 0 0 
= —1/4. 
1 1 
0 0 > 9 
—2t t 
0 0 9 9 


Similarly, the determinant 


1 pti cee (u + 1)" 


1 opti ow. @+ny™ 


can be converted into the determinant of the real solutions we wish to use by multiplying 
by a non-zero determinant. Thus the lemma is established. 
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CorouuaryY. For the solutions obtained in Lemma A, we can suppose that the first 


determinant is not zero. 
Lemma B. Under the hypothesis of Lemma A, every solution 
;x(O), z(1), --°} 
of (6B’) ts a linear combination of the m linearly independent solutions obtained in Lemma 
A. 
Proor. From (6B’) we can infer that there are constants c; such that for p greater 


than or equal to m 


m 


xp) = p c;x(p — 3) 
1 


1 


since the coefficient of A”z is 1. We consider the matrix, 


) 


| x,(0) z,(1) 

( 
! a,(O) z,,(1) 
| | 
| x(0) x(}) 


The above result shows that the (m + 1)-th column of this matrix is a linear combination 
of the first m columns. The (m + 2)-th column is a linear combination of 2, 3, --- , 
m + 1 columns, but since the (m + )-th column is a linear combination of the columns 
1 to m it follows that the (m + 2)-th column is also. The above result can be used to 
show inductively that every column is a linear combination of the first m columns and 
this implies that the rank of the matrix is less than or equal to m. However, the pre- 
ceding corollary now permits us to infer that the rank is m. 

For p greater than m, let us consider the finite matrix obtained by ignoring the col- 
umns for which ¢ is greater than p. This finite matrix is of rank m. From the fact that 
the determinant of order m in the upper left-hand column is not zero, we can infer that 
the last row is a linear combination of the first m rows with coefficients which are de- 
termined solely by the first m columns, i.e. we have 


z(t) = =z c;x,(t) (12) 


for t = 0, 1, --- , p. But since the first m of these equations are adequate to determine 
c; the latter do not depend on p. Since p can be arbitrarily large, (12) holds for every p. 
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SOLUTION OF THE DIFFERENTIAL EQUATIONS OF MOTION OF A 
PROJECTILE IN A MEDIUM OF QUASI-NEWTONIAN RESISTANCE* 


BY 
HARRY POLACHEK 
Naval Ordnance Laboratory, White Oak, Md. 


1. Introduction. An analytic solution for the basic differential equations of exterior 
ballistics has never been found. These equations express the path of a projectile in a 
retarding medium such as air. In practice these equations are solved by a method of 
successive approximations for each new set of constants or initial conditions. This is a 
very laborious process. 

The principal difficulty is the complexity of the law of air resistance. Several solu- 
tions have been found on the basis of simple retardation laws. These, however, are not 
the true laws of resistance encountered by a projectile, as shown by experiment. The 
only solution which is of any practical value to ballisticians is that found by John 
Bernoulli for the Newtonian law of air resistance, R = KV’. This is used in many cases 
where there are no accurate experimental data available. 

Although resistance encountered by a projectile never behaves exactly according to 
Newton’s law, there are many instances where the deviation from this law is “small’’. 
This is true for projectiles moving with velocities less than that of sound. In this paper 
a solution is obtained for the differential equations of motion of a body in a medium in 
which the resistance law is almost Newtonian, or quasi-Newtonian. This is done by 
making use of the theory of differential variations which takes into account differentials 
of the first order only, neglecting differentials of higher order. 

In an illustrative example in the text the author works out a trajectory of a shell 
with an initial velocity of 300 m/s (the velocity of sound is approximately 340 m/s) 
and an angle of elevation of 60°, by both the new method and the standard method of 
successive approximations. The ranges obtained by the two methods are identical to 
the nearest meter, each giving a range of 4240 meters. Since the initial velocity in the 
illustrative problem is almost that of sound and the angle of elevation is rather large it 
may be expected that this method will give accurate results for all trajectories with 
initial velocities below sound. It should be pointed out that the variation of density 
with altitude is also taken into account. 

The derivation of explicit expressions is confined in this paper to the Newtonian law 
of resistance because of its practical usefulness. However, it is apparent from subsequent 
discussions that the methods used are quite general, and are applicable with slight modi- 
fications to any resistance law for which a solution to the equations of motion can be 
obtained in finite form. 

2. Resistance of a moving projectile. If a body is moving in a medium other than a 
vacuum, it is retarded according to some complex law which is related to the density 
of the medium, the size and shape of the moving body, and its velocity. An exact law 
which governs the amount of retardation, or deceleration has not yet been found. The 
motion is influenced by a number of factors which are too difficult for a complete 
theoretical analysis. For instance, when a projectile moving in air exceeds the velocity 


*Received Feb. 14, 1948. 
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of sound, it produces a pattern of shock waves which use up its energy and greatly 
reduce its velocity. Experimental curves are now available which give the retardation 
for several shapes of projectiles and for varying velocities. From these curves it is ap- 
parent that the velocity of sound acts as a dividing line, at which the curve invariably 
experiences a steep and almost discontinuous rise due to the shock waves propagated 
by the moving body. The subsonic region is a region of quasi-Newtonian resistance to 
which the methods of this paper are applicable. In this region the law of resistance 
governing the moving body deviates from Newton’s quadratic law of resistance, R = 
KV’, by a “small” amount. 

3. The equations of motion. For simplicity we will consider a projectile moving in a 
vertical plane, acted upon by the force of gravity mg (where m is the mass of the pro- 
jectile) in a downward direction, and by a retardation force mF in a direction opposite 
to that of its motion. Let us fix a set of axes in this plane, y being the vertical axis and 
x the horizontal axis. The differential equations of motion may then be written in the 
form 


d’x _ pe 
dt ~ dt 
(1) 
dy _pdy _ 
de 84a 9 


where E = R/V. In the general ballistic problem £ is a function of y, the altitude, as 
well as of V, the velocity. EZ is usually written equal to KG(V)H(y), where K is a con- 
stant characteristic of the projectile, G(V) is a function of the velocity alone, and H(y) 
is the density of the air aloft and is a function of y alone. 

4. Solutions of Bernoulli and D’Alembert.’ In 1719 John Bernoulli had obtained the 
solution of the above set of differential equations on the simplifying assumption that 
R = KV”. Later, in 1744 d’Alembert integrated the same equations for the slightly 
more general resistance law R = a + KV”. We are particularly interested in the special 
case R = KV’, but since it does not present any additional difficulties we shall repeat 
here the derivation of the more general solution due to d’Alembert. In Cartesian co- 
ordinates dx/dt = V cos @ and dy/dt = V sin @, where @ is the angle between the tangent 
at any point of the trajectory and the z-axis. Equations (1) may then be written 


< (V cos 0) = —(a + KV") cos 0 

(2) 
a 
di (V sin 6) = —(a+ KV") sin 0 — g. 


Differentiating the left-hand side of Eqs. (2), then multiplying by sin @ and cos @, 
respectively, and subtracting, we obtain 


dt _V 


d@ += g cos 0" 


(3) 


1K, J. Cranz and K. Becker, Exterior ballistics (Vol. I of Handbook of Ballistics), His Majesty’s 
Stationery Off., London, 1921, Ch. IV. 
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It follows immediately that 





dx dx dt Vv’ 

ae ela (4) 
dy _dydt_ _V",. 

do dtde- g “an ¢ 5) 
ds ds dt Vy’ 

do dtd@ gcosé@’ (6) 


Substituting Eq. (3) into Eq. (1) and expanding the left-hand side, we also obtain the 
equation 


g cos 0(V~""")dV — (a+ gsin 0)V “dé = Kdoé. (7) 


This equation may readily be integrated, giving an expression for V in terms of 8, 


V~" = cos” @ exp | ~no" / a(cos @)~* ao 


“ K exp | no | a(cos 6)~* ao 


x i =f a n+1 d6 +> C 
g cos" 6 


For the case a = 0, n = 2 and the initial conditions V = V> and 6 = ¢ we obtain 


(8) 








igs ee ea 
sila 2K cos eC — f(8)] (9) 
1 : do © 
(=a | ara Tor (10) 
I . dé 
° ~ 2K | cos’ 6[C — f(0)] (11) 
l tan 6 d0 
4 = ~ 9K I cos’ o[C — f(6)] (12) 
— C — f(0) 
8 = 2K log Cv — f(¢) (13) 
where 
— teeny 
C > ali 2K(V, cos o)° a S¢) 
and 


f(0) = | sec® 6d0 = }{tan @sec 6 + log (tan 6 + sec 8)]. 


5. Variational equations.” The theory of differential variations was developed during 
the first world war by G. A. Bliss, F. R. Moulton and T. H. Gronwall for the purpose 


*F. R. Moulton, New methods in exterior ballistics, University of Chicago Press, 1926, Ch. IV. 








278 HARRY POLACHEK [Vol. VII, No. 3 


of finding corrections to a trajectory due to non-standard ballistic conditions. A set of 
linear differential equations was derived which had to be solved in order to obtain these 
corrections. This was carried out by numerical methods, in a manner similar to that 
used in the solution of the equations of motion. 
For a given set of initial conditions 
x(0) = 0, y(0) = 0, 
(14) 
2'(0) = 2/, y'(0) = %, 
Eqs. (1) have a definite solution. In terms of ¢ as independent variable, let this solution 
be represented by the functions 
x = $(b), y = V(b), 
(15) 
zg’ = }"(j), y’ = W'(d). 
We now introduce two types of perturbations in the original equations (1): we introduce 
small perturbations in the normal forces acting upon the projectile, thus changing 
Eqs. (1) to the form 


dx a 
;= -E + X 
dt dt 
(16) 
dy , ay , 
S= -E—+Y- 
dt dt g 


where X and Y are the components of acceleration due to these forces, in the x and y 
directions, respectively; at the same time we change the initial conditions given in (14) 
by the addition of small increments a and £ to the values of xj and yf , respectively. 
The new initial conditions thus become 


(0) = 0, y(0) = 0, 


r’(0) = a +a, y'(0) = yo + B. 
The solution to the new set of differential equations (16) with initial conditions (17) 


may be written in the form 


a (ft) + &1), y = V(t) + (0), 
18) 
®’(t) + t’ t), y = W’(t) + 7 (0, 
where é(¢) and y(t) are small compared with ®(t) and W(t). 
Let us now make the assumption that R/V is a function of x’ and y’ but not 
P ly 4: via : ‘ . ’ -f . , i. 
of y; and that it may be expanded as a power series in terms £ = 2 ® and 7 y 


W’. Substituting (18) into (1) and neglecting all differentials of higher order than the 


first, we obtain 


d°® dé OF dk, . 
i + = —| £o 44 rE’ + = 0’ ((P’ + 2’) + X 
ar” 6a 7 ’; ia oY 7 re) + 
(19) 
dv . dn , OE, ., dE, ; . 
+ “1 fet 37 t + The +9)+ 7 -¢@ 


7,2 ( a a = / 
at dt Ox Oy 
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/ 


where E, is evaluated at &’, ¥’. Since (15) is a solution of Eqs. (1) it follows that 





d’& 
—= = —E,®’, 
dt 
| > 
ay a 
a —H,W gq. 


Hence, substituting into (19) and neglecting quadratic terms in ~’, 7’ we obtain 


d’t 4) OM Ob, - 
3 = —Et’ - w| Fee 4 = +x 


dt ~ ? Ox ay’ " 
d’n — ee Ok, e dE, / y 
dr Bon ¥ | 2 dees oy’ ” + 
| or, 
d’t dk, OE, r 
7 = —\E£, ~’ —; t’ — & —>F 7’ X 
dt ( + dx /° oy a+ 
d’n OF, ( : aks) . 
= —wW’ —7 2’ -\E, v’ —> }y’ v. 
rT | ar’ ® + ay)" - 
Setting EK = KG(V) (i.e. H(y) 1), we get 
dk, _ 0G ’ E, ,, 0G 


ar’ ~ Savv =~ VG? av’ 


re) Oe 0G wv’ E, P 0G 
7;= Kk —— = —> ¥ 
ay av V~ Va~ al 


Finally, substituting (23) into (22), 


dt —_ dt’ se Pet +P , + Y 
-:)hlUeS:hCU sid ‘ 


dn = dn’ ™ Q,e’ + Qun’ fe Y 


dt dt 
where 
. , , Ey a 
p, = -(#. + yo =), 
Ek, 0G 


P= = ©! veer’ 


BE, 2) 


Q = —(Bo + Y VG av 


This set of equations may also be written in the form 
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(20) 


(21) 


(22) 


(23) 
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dg =O04+f +040 

dt 

dé’ , . 

dt =0+ P.é +O + Pn’ + X 

(25) 

l 

—1 =0+0 +0+ 7’ 

dt 

dn’ , , ’ i. 

-_—= 0+ Qt’ + 0+ Qon’ + Y. 


The quantities Ey , 6’, G/dV etc. in the expressions for P, , Pz , Q: , Q2 are defined 
for the unperturbed trajectory, which we shall assume has been solved. Consequently 
Eqs. (24) are a set of linear differential equations, the solution to which will give the 
quantities £, &’, », ’ which represent the increments in z, 2’, y, y’ due to a change a, B 
in the initial conditions and the perturbation forces X, Y. 

6. Solution of equations for perturbed initial conditions. We shall first obtain a 
solution to Eqs. (24) for the case X = Y = 0. Under these conditions Eqs. (24) are 
reduced to a system of homogenous linear equations. 

Let us assume that we have found two sets of solutions: 


é’(t) = &(d) 
n’(t) = ni(d 


for the initial conditions 
E’(0) = £(0) 
n’(0) = ni(0), 
and 


e’(t) = £(d) 


n’(t) = (2) 
for the initial conditions 
(0) = &(0) 


n’(0) = 73(0). 


It follows from the theory of linear differential equations that all solutions of Eqs. (24) 
may then be written in the form 


E(t) = vik) + Crkx(0), 


(26) 
n/(t) = Cini() + Con2(d), 
where C, , C2 are constants, provided the determinant 
Hi) CO) | 
D(t) = ~ 0. (27) 


| mi(t) n2(t) | 
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It is also known that the value of this determinant is given by 
t 
D(t) = D(0) exp | / (P, + Q,) at}. (28) 
0 


For the quadratic law of air resistance FE = KG(V) = KV. Equations (25) then re- 


duce to 








df oe"), Kev , 
(29) 
dot _ KW _ er 4. ¥)y 
dt i. KV\1 + a)’. 
These may also be written in the form 
¢/ 
< = —KV(1 + cos? 0) — KV sin 6 cos 61’ 
(30) 


“ = —KVsin @ cos 6’ — KV(1 + sin’ 6)n’ 


where 





__ te\” 1 
V= (z) cos 6[2C — 2f(@))'” (Eq. 9). 


In order to solve Eqs. (30) it is necessary to find a set of two particular solutions to 
this system. Such a set of solutions is in general furnished by the derivatives of the solu- 
tion to Eqs. (1) with respect to the parameters introduced by the initial values. Let the 
solution to Eqs. (1) be given by 


x’ = S(t, x, yd), y’ = W'(t, 2, Yo); 


then 
I’ ov’ de’ ov 
¢’ _ — ; n’ = and t’ _ yy n! as 4 
OXo OXo dYo PYo 


satisfy Eqs. (30). Furthermore, since ¢ is not contained explicitly in Eqs. (1), it may 
be easily verified that 

de’ A 

— rere 








y= 


is also a solution of Eqs. (30). For convenience, we choose the latter solution and also 
the derivatives of ©’, VW’ with respect to the parameter C. These are given in (31) and 
(32) below: 

H(t) = 

, cos 6[2C — 2f(6)]’ 





(31) 
tan @ 





+ 1, 
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with the initial conditions 


K Vo cos $ 








¢/(0) = 
g 
KV; si 
ni(0) = 2 SBE + 3, 
q 
and 
ti dé ] 
E(t) = E(t) a rr 
; Mt ‘ cos 6[2C — 2f(6)]°”~ (2c — 2f(@)]*”* 
(32) 
ie dé _ tan 6 


, 


‘(?) = ;( 2) ——— ———_ — = ae 
™ - I cos’ O[2C — 2f(@)"* [2 — 2f(@)*” 


with the initial conditions 


£3(0) = -(4) ( V, cos ¢)’, 
g 


t— ; 
nx(0) = -(£) (V, cos d)° tan @. 
g i 


We will now evaluate the determinant D(t). We note that 


KV; cos @ — (£) ( V,, cos $)° 
¢ g r\ 3/2 
D(0) = = () (V, cos ¢)”. 
-\ 3/2 g 
Sree ee ; ' P 
KVojsn¢g + 1 - (4) (V, cos ¢)° tan @ 
g g | 


On the other hand 
ie pati Mien Co ae sect 
exp | (P, + Q2) at| = exp (—sx | J dt) = (x) (V, cos ¢) “[2C — 2f(@)]~”. 


Hence, 


D(t) = D(O) exp | | (P, + Q.) ar| = [2C — 2f(0)]"*”. 


0 


It is apparent that D(t) is not equal to zero, except when ¢ = 90°. This case will be 





discussed separately in Sec. 12. 
7. Solution of equations for perturbed forces. We are now ready to solve the set of 
non-homogeneous equations (24). We shall assume that the solution has the same form 
as (26) with the exception that C, , C. are now functions of ¢t rather than constants; | 
namely, 
E(t) = CDE) + Co(OE(O), 
(33) 


Ciba) + Cdn. 


II 


n’(t) 
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Differentiating (33), and remembering that (33) with C, , C, constant is a solution to 
(24) when X = Y = O, we obtain 





rp) Wea 4 gap) MCo _ 
H() + a) 2 = X, 
(34) 
1p) a 4 ry) WC2 _ iy 
n;(t) dt + n3(t) dt =  £ 
Hence, 
, X £3(t) 
a = DO | = f(t) 
| } n3(t) 
and (35) 
| &(t) X 
w= Dt) | = f,(t). 
; | nit) Y 
Integrating we obtain 
Ct) = By + | ft dt =B, + FO) 
(36) 


sf 


B,+ | f(t) dt = B, + Fd) 


0 


C(t) 


II 


where B, , B, are constants which may be determined from the initial conditions. 


Substituting these values into (33) we finally obtain a solution to Eqs. (24), 


E(t) = F\(dE OD + FADED) + Bie) + Beso, 


(37) 
n'(t) = Fi(t)ni(t) + Fld ni) + Bini() + Bond. 
For normal initial conditions B, = B, = 0, hence 
E(t) = FOE) + FeO, 
(38) 


n’(t) = Fi(Oni() + F.(dni(d. 


8. Quasi-Neutonian resistance. Equations (1) express the path of a projectile under 
the influence of a law of air resistance E = R/V. If the resistance is slightly perturbed, 
so that only differentials of the first order must be considered, while differentials of the 
second or higher order may be neglected, the resistance is then given by E + AE. 
When this expression is substituted into Eqs. (1), Eqs. (16) are obtained with 


X = —AEV cos @, 
Y = —AEV sin 6. 


It follows that for the case of quasi-Newtonian resistance Eq. (35) may be written in 


the form 
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cos @ ¢, | 
fi) = —AE V[2c — 2f(@)]*”, 
sin @ = | 
(39) 
| 4 cos 6 | 
f(t) = —AE V[2C — 2f(@)]*”’. 
ni sin 6 | 
Consequently, 
: n§ AE | cos 6 £2 | [2C =e 2f(0@)]'? 
r= | K |. | cost @ OM 
“ | sin 6 n2 | 
(40) 
aml — 88) roe _ on gy 
ie | K ; cos’ 6 o. 
, ni sin 6 | 
We now may write expressions for the increments | 
es [ ss (Ft + Ff) do | 
= Jy (gK)'” cos* o[2C — 2f(0)]'” | 
lt ie an | 
’ = : = FE + F £3 | 
(41) 
see fr ___ (Fini + Fam) do 
1 ~ J, (gK)'” cos® 6[2C — 2f(6))” 
, d / UJ / 
7 =o) = Fant + Font . 


Equations (10), (11), (12), and (41) thus give the complete solution for a trajectory 
with quasi-Newtonian resistance. The procedure for carrying out a solution consists of 
several steps. First we solve Eqs. (10), (11), and (12), using a value of K that would 
give as closely as possible a true estimate of the resistance function. After this has been 
completed the function AE = E — KV is computed. (EZ = KG(V)H(y) is the true 
resistance function, and may in practice also include the variation in density due to | 
altitude.) Then the corrections é, ¢’, 7, 7’ due to a change AE in the resistance function | 
are computed from Eqs. (41). The true values of t, z, x’, y, y’ which we shall denote 
by ¢t, Z, Z’, 7, y’ are then given by the relations | 


t = 
Z=xt+é 
z=2' +?’ (42) 
Y¥-yt 


‘= y' +7’. 


1 
I 
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In order to test the applicability of this method to practical ballistic problems a 
portion of a trajectory based on the standard Army resistance function G, was cal- 
culated by the usual method of successive approximations. The initial values were 
taken as follows: 

Initial velocity = 300 m/s, 

Ballistic coefficient = 1, 

Angle of elevation = 60°. 
The solution obtained is compared in the table below with the solution by the method 
proposed in this paper. There seems to be perfect agreement between the two results, 
as far as the computation was carried. 








Successive Approximations New Method 
t z y z y 
sec. meters meters meters meters 
0.00 0.0 0.0 0.0 0.0 
1.00 146.8 249.5 146.8 249.5 
2.00 287.9 479.6 287.8 479.5 
3.00 424.0 692.0 424.0 691.9 
4.00 555.8 888.0 555.8 887.9 
5.00 683.8 1068.5 683.7 1068.4 
6.00 808.2 1234.5 808.1 1234.4 
7.00 929.6 1386.7 929.5 1386.6 
8.00 1048.1 1525.7 1048.0 1525.5 








9. The adjoint system of equations. In many ballistic problems only the terminal 
conditions of the trajectory, such as the range of the projectile when it strikes the 
ground, are required. In that event it is not necessary to obtain the path of the complete 
trajectory. By solving the system of equations adjoint to (25) it is possible to obtain 
directly the terminal values and thus to eliminate a large percentage of the computa- 
tional work. In addition to this gain the method of the adjoint system yields expressions 
that can be used to find the effects upon the normal trajectory due to non-standard 
conditions of the air, such as wind or temperature variations. There is a third advantage 
to the computer in this method in that the final results are given for equal intervals of 
6, which is not true in the case of the formulas derived in the preceding section. 


The set of equations 


_ dp 


di 0+ 0 +0+0 


II 





ap p+ Pip’ +0 + Qe’ 
(43) 


“ 0+0 +0+0 


do’ 


dt 


0 + Pip’ + o + Q,0’ 
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whose coefficients are the same as those of Eqs. (25) with rows and columns interchanged, 

is known as the system of equations adjoint to (25). If a solution to (43) is found for 

the variables p, p’, a, o’ the quantities é, ’, n, ’ may then be expressed in terms of these 

variables. The relationship between these two sets of variables may be easily derived. 
Multiplying Eqs. (25) by p, p’, o, o’ respectively, and adding we obtain 


dé ,d , dy 


: - , dn’ 
Pat t p dt . dt 


5 is dt 


(44) 
= §’(p + Pip’ + Qio’) + n’(o + Pop’ + Qro’) + (p’X + o’Y). 


Substituting Eqs. (43) we obtain 


lé lt’ d ln’ 1 lp’ do’ y , 
ag , UB 7 , GN ge t. , Up 4 n+ n’ lo’ _ p’X + o’} (45) 


ctr 


eee ile” ei ee dt * "¢ dt 
or 
d & sy / / / 4 , y » 
7, PE + P'S + on + o’n’) = p'X + 0'Y. (46) 
Finally, 
T oT 
| 0 + pt’ + on + a | sas | (p’X + oY) dt (47) 
0 “Q 


where 7’ is the time at the terminal point of the unperturbed or Newtonian trajectory. 

10. Solution of the Adjoint System. It follows immediately from the first and third 
of Eqs. (43) that p = constant and that o = constant. We shall denote these constants 
by p and ag, respectively. The second and fourth equations may then be written 


dp’ os , 
_ = p+ Pip’ + Qi 
dt 
(48) 
do’ ae a P 
— — =a Pr.’ + Q.0’, 
dt 
or 
dp’ ( : » By eG E, 0G 
- =p-—-Ii, eo” — -|p’ — P’’ — + 
aa” + ® Vaav/? VG av’ 
(49) 
da’ E, oG ( 7 * 20) 
-- : — P/V’ — -p’ —\E, yy’ —— —— jg’. 
dt ’ J G ol F . } G ot “ 
It may be verified that one solution of this system is 
pi?’ +4,’ + pio’ + of” = —-K, (constant). (50) 


This follows from the fact that § = %’, = ©”, 7 = W’, 7’ = W” is a solution to Eqs. 
(24), and from (46) which reduces to 


d 


xT ( pé ot pt’ ot on te on’) _ 0 for t => Q, Y — 0. 
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Transposing the terms of Eq. (50) and dividing by ®” we obtain 


— ar — wy ayer 
—pi _ piP’ + oe oiV + K, (51) 





which may be used with the second equation of (49) to eliminate p{ ; thus 


Ey 9G p.® +a, + of ¥"” + K, 
VG av p”’ 


1, , 
—~- 5,40’ 
dt 





(52) 
» E, 0G 
—\E vv’ ——le . 
( ot VG av/ 
Remembering that ®’ = 2’, VW’ = y’ of the undisturbed or Newtonian trajectory, and 
that 6” = —£,%’, the preceding equation reduces to a linear differential equation of 
the first order which may be solved immediately; namely 


Wi + A(oi + BY) = 0 (53) 


where 


gv’ ome .. E, , 


A(t) VdV 


II 


d log G 


VavV (p,%’ + o,¥’ + K,). 


B(t) = G} ~—- Vv’ 


The solution of this equation is 


o{ = exp (- [ A(t) an) c. = [ B(t) exp ( [ A(t) at) at}. (54) 


“0 “0 


The expression for p{ is given by Eq. (51) above. 
We shall want only the solution for the case of = 0 at t = T. Hence Eq. (54) may 
be written 


at 


o(f) = exp (— | 
| [ B(t) exp ( [ A(t) at) dt — [ B(t) exp ( [ A(t) at) at}. 


0 0 


A(t) at) 
(55) 


We are particularly interested in the solution of the above equation for the special 
case FE = KG(V) = KV; hence d log G/VdV = 1/V”’. From Eas. (3), (4), (5) and (9) 
it follows that, 


A(t) _ gsin @ — KV 


j 
and 
" sna en ~~ | cme KV --f = a 
| A@ at = / coe t= — | tan odo + a sa aC TOI ° 
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This may be integrated, giving the solution 


Ad) dt = low 128-2 | C= SO |". 
[ A(t) dt = log {see E _ fo | h 


/0 


and finally, 


al f° ssa el ome | g= fe |” _ («) ——— 
— (/ A(t) it) = cosé LO — f(a)1 ~\K) Voces o2C —2fo)7° 

















Also, 
ee neg . 
Bit) =o, -— — V (p,V cos 6+ o,V sin 0+ K,) 
= 5, cos 6 — pf, cos Osin 0 — ins j 
Hence, 
; x 1 oe i tm K, sin 0 
pty ono (A092) a= gr tor f° [ (Ep utin 
| ()exp\ {, AW dt) dt = er ore i g/ PC — 27)” 
_ % cos 0 = 20s 2] 
20 — 2f(6) 6. 
Finally, 
= hay BO = shor” f° is _Kisin@ __ @ cos 8 ord do 
1 (Kg)'” cos r) Je q [2c A 2f(0)]'”7 IC rae 2f(0) 
(57) 
,_ _ PO’ +a’ + oi ¥"’ + Ki 
a= @” 


where w is the angle at the terminal point of the trajectory. 

11. Range and time of flight. As is well-known, the solution derived for the adjoint 
system of differential equations (43) may be used to obtain the range, the total time of 
flight or other terminal values of the perturbed trajectory. The quantities that are of 
greatest interest are the range (the horizontal distance from the muzzle of the gun to 
the point where the projectile strikes the ground) and the time of flight. 

Let us suppose that Eqs. (10), (11), and (12) have been solved for a trajectory 
with given initial conditions. Let 7 be the time of flight to the ground, and X,7 the 
range. These equations are based on the Newtonian resistance law, and the resulting 
trajectory is an approximation to the true trajectory. Since the true trajectory which 
obeys a quasi-Newtonian resistance law is not much different from the Newtonian 
trajectory we shall consider only differentials of the first order. In Fig. 1, Xz is the range 
of the Newtonian trajectory, X7- is the range of the quasi-Newtonian trajectory, P 
is the location of the projectile at time 7, and PO is a line normal to the z-axis. The 
change in horizontal range ¢(7) = X70; while X70 = —2’(T)n(T)/y'(T). Hence 
X7Xyr- the total increment in range, is given by the quantity (7) — 2’(T)n(T)/y'(T). 
Likewise, the increment in time of flight, which is the time the projectile moves from 
P to Xr. , is given by —n(T)/y'(T). 
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In the preceding section we have found the solution to the adjoint system of differ- 
ential equations. We may now assign initial conditions to p, p’, o, o’. Let us first assign 
the set of initial conditions 
«'(T’) 
y(T)’ 


Substituting these in Eq. (47) we obtain, since (0) = (0) = 0 


AT)=1, p(T)=0, off) = o(T) = 0. (58) 


os 2) oes Pe a 
ECD) — Fy WT) = w'OR'O) + 0' On’) + | (o'Xitie’Y) dt. (69) 


QUASI-NEWTONIAN TRAJECTORY 






NEWTONIAN TRAJECTORY 





Fia. 1 


It will be noticed that the left-hand side of Eq. (59) is the increment in range de- 
rived from geometrical considerations. Thus Eq. (59) gives an expression for the difference 
in range between Newtonian and quasi-Newtonian resistance in terms of the functions 
o’ and p’ of the adjoint system. These are given by Eqs. (57) above. 

In order to obtain the increment in time of flight we assign the initial conditions 


1 





AT)=0, p(T) =090, . o(T) = Wm’ ° (T) = 0. (60) 
These give, by substitution in (47) 
(T Rn ere i ita ” pe laa 
— Teg = POO) + o'O)n'O) + | (o'X + 0°Y) dt, (61) 


which is an expression for the increment in time. The functions p’ and o’ may again 
be obtained from Eqs. (57). It will be noted that the value of K, is not the same for 
range and time of flight but must be determined from Eq. (50) before solving Eqs. (57). 
For the increment in range K, is equal to zero, while for the increment in time of flight 
K, = 1. 

A complete trajectory with the initial conditions used for the illustrative problem 
in Sec. 8 was computed by the standard method and by the method of the adjoint 
system. The ranges obtained by the two methods were identical to the nearest meter, 
each giving a range of 4240 meters. The actual values obtained were 4239.8 meters and 
4240.2 meters. This may be considered complete agreement since the inaccuracy due 
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to the method of successive approximations is of the same magnitude as the discrepancy. 
For practical ballistic problems this is far greater accuracy than is usually desired. 

12. Vertical trajectories. It was pointed out in Sec. 6 that the equations derived 
there for a trajectory with perturbed initial conditions do not hold for a vertical tra- 
jectory, since for such a trajectory dx/dt = 0. In this section we will derive simplified 
formulas for this special case. It will be of interest to note here that for the case of a 
vertical trajectory the solution of the equations of motion for Newtonian resistance 
with variation in air resistance due to an exponential law of density aloft is known. 





The solution may be given in terms of the function Hi(z) = f2. (e*/z)dz. 
In one dimension the equation of motion is 


dy _ py 


= , (62 
dt dt g ” 


Let the solution to this equation for the initial conditions 


y(0) = 0 


y'(0) = 


be given by 


Also, let the solution for non-standard initial conditions 


y(O) = 0 


I 


y'(0) =y+8 


be given by 


y = V(b) + nb) 


II 


y’ = V(t) + n'(0). 


If we assume in addition that EK = R/V is a function of y’ = V alone and may be ex- 
panded in terms of n’ = (y’ — W’), we may write 
d’y . dn | : dE, ‘ " 
~ale —, = —| £, "10 — g. (63) 
dt’ + dt’ + dy’ | +- | ) g 


Substituting Eq. (62) into (63) and neglecting terms in 7’ we obtain 


dn dE, ; 
Sm ~ _— ! , ! e (64) 
dt” ’ y dy ’ 
Since y’ = V and E can be written E = KG(V), the above equation reduces to 


d'n _ dy _ -| sakes me ag | : en 
a 7. = KG(V) + KV dv |"° (65) 
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Hence 


f=Cem | -/ (Kav) +KV ey ar}. 


For a small abnormal force Y acting on the projectile, Eq. (62) may be written 


dy_ _,d _ 
a fat -& 


For this case . 
d’n _ dn! _ | < 7 a ’ 4 
de” ae = [RAV + RV GG WW + F 


which is a linear differential equation. The solution is given by 


n’ = exp | -x | (av) +V a at| 


Ac, 4 / Y exp | x / (av) 4 y 2) at| at 


For the special case G(V) = V 


f= ep (—2Ky)| C, 4 | Y exp (2Ky) at}. 
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(66) 


(67) 


(68) 


(69) 


(70) 
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SOLUTION OF STEADY STATE TEMPERATURE PROBLEMS WITH THE 
AID OF A GENERALIZED FOURIER CONVOLUTION*' 


BY 
A. W. JACOBSON 
Wayne University 


1. Introduction. The purpose of this paper is to illustrate the use of the finite Fourier 
transformation and a generalized convolution in the solution of steady state boundary 
value problems. Solutions of some basic problems are given in terms of special functions 
introduced, and by means of a modified Duhamel formula other solutions in turn are 
expressed as functions of the basic solutions. 

The finite sine and cosine transformation are defined by 


Il 
_ 
~ 
bo 
~ 
~— 
~ 


s{r(a)} = f ‘Adtemebonf{M & 


C{F(az)} = 4 F(x) cos nz dx = f,({n) (n = 0, 1, 2, ---), 


respectively. These transformations applied to the derivatives of F(x) yield, for example, 
the following: S{F’(xz)} = —nC{F(x)}, S{F’(z)} = —n’S{F(z)} + n[FO) —- 
(—1)"F(x)], and C{F’(xz)} = nS{F(x)} — FO) + (—1)"F (x). 

If F(x) in (—2z, 2x) and G(z) in (—z, 7), are sectionally continuous functions, then 


the function 


F(z)*G(2) = [ F@ - WOW) dy (1) 


is called the convolution of F and G in the interval (—z, 7). 

If F(x) is an odd and G(z) is an even sectionally continuous function and if 
F(x + 2x) = F(z), then the product of the transforms can be written in terms of the 
transform of the convolution, for example,’ 

S{F(x)}C{G(z)} = 2 S{F(x)*G(a)}. (2) 

2. A generalized Fourier convolution. Let F(z, y) be a sectionally continuous 
function of z and y in the square 0 < x < 7,0 < y < =z. A generalized convolution 
F*(x) of F(x, y) corresponding to the iterated finite sine transformation 


S{S{F(a, y)}} 


| [ F(x, y) sin nx sin n'y dx dy 
“0 “0 


f(n, n’) (n,n’ = 1, 2, --:) 


is defined by 


er 


F*(z) = — | F(x — y, y) dy, 





*Received August 24, 1948. 

1Presented to the American Mathematical Society April 15, 1948. Also see abstract in Bull. Am. 
Math. Society, vol. 54. The author expresses his thanks to Professor R. V. Churchill for many valuable 
suggestions. 

2R. V. Churchill, Modern Operational Mathematics in Engineering, McGraw-Hill, 1944, p. 274-276. 
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where F(z, y) is an odd periodic extension of F with respect to x and an odd extension 
with respect to y. 
When n’ = 2, it can be shown’ that 
S{S{F(z, y)}} = 2 C{F*(a)}. 
In case F(x, y) = F(x)G,(y), the function F*(z) is the convolution (1) of F, and G, . 
3. Special functions. In this section two key functions are introduced as well as their 
sine and cosine transforms. In the sequel solutions of some basic boundary value prob- 


lems are expressed in terms of these functions. 


If |r| < 1, then log (1 + re’) = — >>%_, (—r)’e'”’/». Also if R? = (1 +r cos 6)? + 
r’ sin’6 and tan ® = rsin 6/(1 + r cos @), then 
log (1 + re’’) = log Re** = log R + i@ (5) 
Equating the imaginary parts of (5), we get 
rsin@  _ = (—r)” sin v8 
oes YT he > v ; 


When |r| < 1, the series is a Fourier sine series. The sine transforms of the function 
are, except for the normalizing factor 7/2, the Fourier sine coefficients, i.e., 


rsin x a (—1)’r" 
S\ arctan ————p = —2 A < 
s{arctan 1 +7 cos +} > wn * Ir] <1, 


where 6 has been replaced by x. Since (—1)"*’f,(n) = S{F(x — z)}, then 


rsin x xr” 
SS arctan ———__/ = =-—. 
l1—rcosz 2n 


Setting r = e *, the above formula becomes 


S{arctan 1 ? w= >0. (6) 
e* — coszx 2n 
Let b,(n, u) denote the function 
Sih te = + texp {=[(2v + Ix — un} — exp {-[2 + Ix + ula}. @) 


Let the inverse sine transform of the function bo(n, u), S~*{bo(n, u)}, be Bo, (x, u), then 
according to formula (6), when r — u > 0, 





° sin x 
Bo(z,u = arctan exp {(2v + l)r — u} — cosz 


— arctan ain & | (8) 
arcté : 
exp {(2» + l)r + u} — cosz 





* t sin z sinh u 
arctan = ; 
cosh (2v + 1)r — cos zx cosh u 


y=0 





A. W. Jacobson, A Generalized Convolution for the Finite Fourier Transformations, Thesis, 
University of Michigan, 1948, p. 9-12. In the sequel reference to this thesis will be marked [3]. 
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Equating the real parts of (5), we have 





- a _@ee(-DT __ 
log 1+ 2reosx+r° . > eee 


Hence the cosine transformation yields 








f \  (=—1)"r" 
Cs —_— a ? eee = ~ * 
log 1 + 2recosx + ef n P . 1,2 
= n=0O0 
Since (—1)"f.(n) = C{F(a — z)}, and setting r = e *, we get when z > 0 
cf tog Sr Rang — >} = = : n= i1,2,--- 
4|COsn Zz — COS TZ n 
) (9) 
= 0 n=O 


Let b,(n, uw) denote the function 


cosh nu — | 
7 — = >> = fexp {—[(2v + 1)r — ujn} + exp {—[(2v + Ix +uln}}. (10) 
n sinh nx —F 
And let B,,(x, u) denote the inverse cosine transform C™'{b,(n, u)} of the function 
b,(n, u). When x — u > 0, according to formula (9), we get 


ee ee il - exp 202 + Ir} 
a | 4[cosh (2vr + r — u) — cosx][cosh (2yr + x + u) — cos az] 


vy=0 





(11) 
4. Problems in two dimensions. Basic problems (A). Let Uo(z, y) be the solution of 
the following steady state temperature problem in the region R:0 < 2 < 7,0 Sy < =: 


7U, ?U, 
al 3 = 0 in i, 


dx Oy” 
U,(+0, y) = Uo(x — 0, y) = 0, 0O<y <r, (A) 
Uz, +0) ==—*,  Ulz,r -0 =0, O<2<r. 

Tv 


The sine transformation of problem (A) with respect to x yields 


duo 


dye 7 tol v) = 0, 


(A’) 


l 
un, 0) = ” U(n, r) = 0. 


The solution of the transformed problem is 


sinh n(x — y) 


U(r, y) = = 
n sinh nr 
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In the previous section, formula (7), this function was designated by the function 
bo(n, u), 1.€., 


Uo(n, ¥) = bo(n, e — y). 


According to formula (8), its inverse sine transform is the function Bo,(z,  — y), so 
that 
U(x, y) = Boi(x,  — y). (12) 
Problem (B) 
aU , #&U , 
az? ay? = 0 in i 
U(-+0, y) = Ue — 0, y) = 0, 0<y<r, ” 
U(x, +0) = H(z), U(z, x — 0) = 0, O<2<-e. 
The transformed problem is 
2 
a — n'u(n, y) = 0, 
y (B’) 


u(n, 0) = h(n), u(n, r) = 0. 


Upon multiplying the equations in problem (A’) by the function nh(n) of the param- 
eter, it is evident that the product nh(n)uo(n, y) is also a solution of problem (B’). If 
the solution of the latter is unique then 


u(n, y) = nh(n)ud(n, y). 
Since U,(0, y) = Uo(x, y) = 0, then nS{U(z, y)} = C{d/dax U,(a, y)}. Hence 
S{U(a, y)} = S{H(z)}nS{U,(z, y)}. 
becomes 


S{U(z, } = simaicy2 Udle, w}. 


The product of the transforms here can be written, according to formula (2), as the 
sine transform of the convolution of the two functions, and making the inverse trans- 


formation, there results 
: Pe R 
U(x, y) = 4H(z)* - U,(a, y). 


Since U,(z, y) = Bo:(z, r — y), equation (12), the last result becomes in terms of the 


convolution integral (1) 


U(z, y) = ; [ H(z — ») = oe ey 


In view of formula (11) 8/8 Bo,(A,  — y) is an even function of A, hence H is to be 
extended to (—7z, 0) as an odd function. 
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Problem (C). 


92 J 93 J . 
aU al = F(z, y) In R 


ox” Oy’ 
U(+0, y) = U(r — 0, y) = 0, O<y<r ©) 
U(z, +0) = U(z, x — 0) = 0, 0O<a2<r. 
The sine transformation of this problem with respect to x yields 
d’u 
—_ n'u(n, y) = f(n, y) 
Y (C’) 


u(n, 0) = u(n, r) = 0. 


The solution of the transformed problem, in terms of the Green’s function g(n, y, »), 








is 
win, ») = | alm, yy ofr, w) di (13) 
where 
sinh nu sinh n(x — 4 
g(n, Y; h) = — ne y) ed < y 
sinh ny sinh n(x — yp) 
n sinh nr : way, 
or 


cosh n(x — y — ») _ cosh n(r — y + ») < 
2n sinh nr 2n sinh nr , esYy 





gn, y, uw) = 


cosh n(r — y — 4) _ coshn(r —u+y) 
2n sinh nr 2n sinh nr 





; uD Yy.- 





According to formula (10) this can be written 


g(n, y, w) = 4b,(n, r — y — w) — 3D(n, er —y +B), assy 


=4bh(n,r-—y-—u—3dm,r-wt+y), w2y. 
The inverse cosine transform of b,(n, u) is, in view of formula (11), the function B,,(z, u). 
Hence 


G(z, y, uw) = 4 BiA(z, 7 —y — w) — 3 Bult, er —y + yw). (14) 


where 
C"{g(n, y, w)} = G(z, y, x). 


Writing solution (13) as 


stu, a} = [ CtGe, v, o)SIF, w)} du, 
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and replacing the product of the transforms by the sine transform of the convolution 
and making the inverse transformation, we obtain 


o* 


‘ l ; . 
U (2, y) 9 | CAL, YB, u)*l(a, uw) du, 


or 


l a , 
U(z, y) = = | | G(x — A, y, mF, w) dd du. 


“0 . ' 


Substituting for G from (14), there results 


U(x, y) = ; | | [Bia(a \, 7 y — p) B(x — A,r — y + wIFA, pw) dr du. 
The function B,.(z, uv) is an even function of x, hence H(z, y) must be extended to 
(—7, 0) as an odd function of x. ! 

The above solutions can be shown to satisfy the respective boundary value problems 
for a wide class of functions involved. See [3] p. 58-75. 

5. Resolution of boundary value problems with the aid of a generalized Fourier 
convolution. Let U(z, y) be the steady state temperature in the region R: 0 < x < fg, 
0<y< 7,0 <2 < a, satisfying the following conditions: 


VU vu au . 
( . = F(z, y, 2) in R, 


Ox Oy . dz” 

UO, y, 2) = H,(y, 2), U(x, y, z) = 0, (D) 
U(z, 0, 2) H.(x, y), Ue. 2,2) = 2, 

U(a, y, 0) H(2x, y), U(x, y, x) = 0. 


Three of the boundary conditions have been written here as homogeneous since this 

involves no loss in generality. The sine transform of this problem with respect to x is 

ou ou , 
+ = f(n, y, 2), 


-n'u(n, y, 2) + nH (y, 2) +3 + 5 
; OY Oz 


u(n, 0, z) = h(n, 2), ulin, x,z) = 0, 


u(n, y, 0) = h;(n, y), u(n, y, r) = O. 


Now let V(2, 2’, y, 2), depending on the parameter x’ independent of x, y, z, be the 


solution of the following auxiliary problem: 


o°V a’V oV T—- 2 


: = F(x’, y, 2) in R, 
Aa oy” t Oz T y 
: rs ‘ - 
V(O, x’, y, 2) = H,(y , 2) Vir, 2’, y,2z) = 90 
(I) 
T , , 
Vix. 2’, 0,2) = H,(z’, 2) Viz, 2’. 8,2) =0 
T 
’ _ toms ; , 
Viz, 2’, y, 0) = Hix’, 4), Viz, 2’, y¥, 7) = 0 
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The sine transform of (E) with respect to z is 


— zg’ 9 a’ l 
no(n, x’, y,z) +n ——=— ily, 2) -; : oy? sk “= =. F(z’, y, 2), 
us OY 
, l , , 
v(n, x’, 0,2) = = H(z’, z), o(n, 2’, x, 2) = 0, 


0. 


II 


] 3 
v(n, xz’, y, 0) = . H(z’, y), v(n, x’, Y, 7) 


Next let #(n, n’, z) be the sine transform of v(n, z’, y, z) with respect to 2’. When n’ = n, 


we get 
—n’o(n, y, 2) + H,ly, 2) + aa ov 3+ a. =- 3 Y, 2), 
v(n, 0, z) = : h,(n, 2), v(n, x,z) = 0, (E’) 


0. 


o(n, y, 0) = : h;(n, y), v(n, y, ™) 


If we now multiply the equations in (E’) by the parameter n, we note that problems 
(D’) and (E’) are equivalent, and 


u(n, y, Z) = nv(n, y, 2). (15) 


Since i(n, y, z) is the iterated sine transform of V(z, 2’, y, z), it follows, according 
formula (4), that 


u(n, y,z) =n ct} V*(z, y, a}. 


- -si2 — V*(z, y, a}. 
And therefore 
} a 
U(z, y,2) = 9 ap VW) 
which can be written (formula 3) 
2 ; 2 
U(z, y, 2) = of Viz — x’, x’, y, 2) dr’. (16) 


Formula (16) is the Duhamel integral extended from time to space coordinates. 
Problem (D) can further be resolved into simpler problems. We first note that 
problem (EF) can be written as 


V(x, 2’, y, 2) = oa Vi(z’, y, 2) + V.(z, 2’, y, 2) + V3(2, 2’, y, 2), 
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where V, , V. and V, are solutions of the problems: 
139 2 3 


av; a’V, - 


= ( 
oy” Oz” ), 
r , , , , (a) 
V,(z’, 0,2) = H(z’, 2), V,(z’, x, 2) = 0, 
V(x’, y, 0) = H(z’, y), V(x’, y, ) = 0; 
a’V, av. a’V. 
— — =~ —_;- = ( 
Ox” + oy” Oz ' 
V.(0, 2’, y, 2) = aes H(y, 2) — Vi(z’, y, z 
= Q(z’, y, 2); (b) 
V(r, 2’, Y; 2) = 0, 
V.(z, 2’, 0,2) = V,(z, 2’, x, 2) = 90, 
V(x, x’, y, 0) = V2(z, x’, y, r) = 0; 
and 
V, , Vs, , OV, _wr-tpH, 
Ox” ay” é27 a P(x’, Ys 2), 
V;(0, x’, y, 2) = V3(x, x’, y,2z) = 9, (c) 


I 


V,(2, 2’, 0, Z) V3(z, x’, us 2) aid 0, 


V3 (2, 2’, Y; 0) _ V3(z, of Y; 1) _ 0. 
It can be readily shown, see [3] p. 19-23, that the solution U(z, y, z) of problem (D) 
then takes the following form: 
—_— Lo vel, , 10 vs . 
U(z, y, 2) = Vi(2, y, 2) — 2 az V3 (x, y, 2) — 9 az V# (za, y, 2). 


The solution V, of the two dimensional problem (a) is the sum of solutions of the 
type of problem (B), Sec. 4, that is, 


V,(z’, y, 2) = : H,(x’, z)* A Bo ly, « — 2) 
2 Oz 


i , 
+ 5 Hele’, )* a Boy,  — 2) 


We shall next proceed to obtain the solution V, of problem (b). Let v.(z, x’, m, 2) 
be the sine transform of V,(z, x’, y, z) with respect to y and also let #,(z, x’, m, p) be 
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the iterated transform of V,(z, x’, y, z) with respect to y and z. Then applying this 
iterated transformation to problem (b) we get 


DD. rT} , , 
oa — (m’ + p’)i.(z, x’, m, p) = 0, 
(b’) 


v(0, x’, m, p) = G(x’, m, p), (7, x’, m, p) = 0, 
where g(x’, m, p) = S{q(z’, m, z)}. 
We shall now formulate the following basic problem (d) in three dimensions, and 
then express the solution V, in terms of its solution: 
a°l Jo a” Uo a” Uo 


+ Sa + a = 0, 


Ox” ay” Oz” 





I 


= 2 Fs = 
r ’ U(r, Y; z) 0, (d) 


T 


U,(0, yy; Z) 


U(x, 0, 2) U(x, x, z) = 0, 


I 
= 


U,(z, Y; 0) _ U,(2, Y; 7) 


Let u(x, m, z) and %(x, m, p), respectively, be the sine transform of U,(z, y, z) 
with respect to y and its iterated transform with respect to y and z. Then 


Oo te? a te - 
ax (m° + p )uo(x, m, p) = 0, 

(d’) 
u(0, m, p) = ide Uio(x, m, p) = 0. 


mp’ 

If we multiply the equations in problem (d’) by the function of the parameters 

mpq(x’, m, p), it is evident in view of problem (b’) that 
v(x, x’, m, p) = G(x’, m, p)mp u(x, m, p). 
That is 
S{v2(x, x’, m, z)} = S{q(2’, m, z)}mpS{uo(z, m, z)}. 

Since uo(x, m, 0) = uo(z, m, r) = 0, then pS{uo(x, m, z)} = C{duo/dz (x, m, z)}. Hence 
the product of the transforms can be replaced by the sine transform of the convolution 
of the two functions. Upon applying the inverse transformation, we get 


l 0 
v.(z, x’, m, 2) 5 q(x’, m, z)*m 3 Uo(x, mM, 2) 


I [ Gy) 
= — g(x’, m, z — A)m — u(x, m, A) Ad, 
a]. q(x’, )m ~~ wol ) 
which in turn can be written formally as follows: 


aT 


ss 1 a -) = 
S{V.(z, z’, y, z)} = =| S{Q(z, y,z — A)}m s\2 U,(a, y, wf dn. 
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Since again 0/0 U,(2, 0, A) = 0/dA U,(x, x, A) = 0, then, repeating the steps above, 
we obtain 


2 


o 7 Fs 
[ [ Q(x, y — wz — ®) HG Vola, m, d) dd du. 


vg du 


| 


Vz, x’, y, 2) = 


The solution U,(z, y, z) can be obtained with the aid of the Fourier sine series. 


The solution of problem (d’) is 
1/2 


sinh (r — x)(m? + p’) 
mp sinh x (m* + p’)'” * 





U(t, M, p) = 


These functions are, except for the constant factor 1/2, the Fourier sine coefficients for 
the function u(x, m, z), which in turn are in like manner the sine coefficients for the 


function U,(z, y, z). Hence 





4x > sinh ( x—2(m+p)'” 


U,(z, y, 2) = waka +a sin pz sin my. 


m=l1 p=l 


Solution of problem (c). The iterated sine transform of (c) with respect to x and y is 
Fie = 1 

<3 — (n® + m\i,(n, 2’, m, 2) = = f(z’, m, 2), 

Ox n (c’) 


v,(n, x’, m, 0) = 0,(n, x’, m, r) = 0. 


In terms of the Green’s function g(n, m, z, 4) the solution of (c’) has the form 








b(n, x’, m, 2) " gn, m, 2, aoe m, p) du i 
where 
sinh p(m? + n’)'” sinh (z — 2)(m? + n’)'” 
gin, m, 2, w) nm sinh a(m 4. ny? ’ BS 2Z, 


and with » and z interchanged when yu > z. 
Performing the iterated inverse sine transformation with the aid of Fourier sine 


series, we get 


> v,(n, x’, m, z) sin nz sin my, 


m=1 


V,(z, 2’, y, 2) = 


“| 
- 


~ 


n= 


where #,(n, x’, m, z) is defined by (17). 
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GENERATION OF SURFACE WAVES BY A MOVING PARTITION* 


BY 
E. H. KENNARD 
David Taylor Model Basin, Washington, D. C. 


Artificial waves in water are frequently generated by imparting rhythmic motion to 
a boundary of some sort. Two cases of wave prodyction in deep water have been treated, 
for harmonic motion only, by Havelock.’ The more interesting of these cases will be 
treated here by a different method, and expressions will be obtained for the initial 
phases as well as for the permanent regime; in addition, the solution will be extended to 
liquids of limited depth. 

The harmonic case requires, for the most part, another two-dimensional solution of 
the Laplace equation, subject to a special set of boundary conditions. Besides the familiar 
condition at the free surface expressing the action of gravity, and perhaps a condition 
on a parallel bottom, the horizontal motion is here prescribed over a vertical boundary. 
The functions found are interesting in that they constitute an infinite variety of char- 
acteristic solutions, all asymptotic at infinity to the same functional form. This asymp- 
totic form is the sine function that is familiar in the treatment of gravity waves and 
represents the only bounded solution when the free surface is unlimited in extent. 


[. Case oF INFINITE DEPTH 

Consider a semi-infinite body of incompressible, frictionless liquid having a free 
surface and of infinite depth. Let it be limited at the left by a boundary that executes 
infinitesimal displacements from an initial vertical plane position. Initially at time 
t = 0, let the liquid be at rest with a horizontal free surface. 

Draw the y-axis downward in the initial plane of the vertical boundary and the 
x-axis perpendicular to this plane and in the horizontal plane that is occupied initially 
by the free surface. Then the motion excited by the motion of the boundary will be 
two-dimensional, occurring in planes parallel to the xry-plane, so that nothing is a 


function of z. 
Let the infinitesimal horizontal displacement s of a point of the boundary at time ¢ be 


s = SF(y, 0), (1) 


where S is a constant having the dimensions of length. Thus the function F(y, 0) is 
dimensionless. Initially, F(y, 0) = 0. It will be assumed that the function F(y, ¢) is 
sufficiently smooth and vanishes sufficiently rapidly as y —© to justify the operations 
that are to be performed. 
Since the liquid is assumed frictionless and incompressible, there will be a velocity 
potential satisfying Laplace’s equation, 
<4 So eG, 


ax ay" 


IV 


A solution is required for x = 0, y = 0, subject to two boundary conditions. Over the 





*Received June 21, 1948. 
1Havelock, Phil. Mag. 8, 569 (1929). 
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plane of the moving boundary the vertical component of velocity of the water is un- 
restricted, but the horizontal component must be the same as that of the boundary. 


Hence for x = 0, 


_§ _% _ 8 _ op 
= — 5 = 5, = SFLy, d- (2) 


At the free surface the boundary condition for infinitesimal motions may be written’ in 
terms of the elevatfon 7 of the surface: 


= 1 (2) on _ (28) 
n= 2(%) (3a), an (#) (3b) 


where g is the acceleration due to gravity. 

A solution of the differential equation satisfying these boundary conditions will be 
constructed in two steps. For the moment, imagine gravity to be absent. Then, since 
the displacement of the boundary remains small, the motion of the liquid can be re- 
garded as due to a suitable distribution of horizontal line sources along the positive 
y-axis, continued above the free surface along negative y as negative image sources in 
order to keep the pressure constant on the free surface, or in the zz-plane. The velocity 
potential due to a line source located at (0, +8) can be written —A log r? = —A log 
[z? — (y +8)’], where A is a constant and the flux toward one side is 27A; the desired 
flux from dy, however, is v,dy = SF,(y, t) dy by (2). It is thus found that the desired 
velocity potential ¢, and the corresponding horizontal velocity component v,, are 

S ‘ia x + (8 + y)” 


= 3 | Fi, 0 log es as, (4) 





_ — % _ 8 [7 py ta ,| 
— id | Fs, | zap 7+e+yl™ 








where F,,(8, t) = dF (8, t)/dt. 
It is readily seen that ¢, is a solution of Laplace’s equation and that d¢,/d¢ = 0 when 


y = 0. Consequently, the pressure remains zero at y = 0, as it must for infinitesimal 
motions in the absence of gravity. As x — 0, the second term in the expression for v,, 
vanishes. In the first term, however, the integrand is small except near y = 8, so that 
in the limit as z — 0, F,(8, t) may be replaced by F,(y, ¢) and the lower limit of the 
integral may be replaced by — ©. Thus 





, S 26 1g 
lim v. = ~ Fly, 0) i 3 Ta ra SF,(y, 0), 


z—0 


as is required by (2). 
The corresponding upward velocity of the surface of the liquid is 


w= (%) 8 ng, 9 22 ; 
vn = (% a “ i F,(8, DP ae: (5) 





The effect of gravity may now be introduced by the following device. During each 
element of time dt the motion already assumed adds an increment to the elevation of 


*Lamb, Hydrodynamics, Cambridge University Press; 6 ed., Sec. 227. 
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the surface represented by »,, dt. In the absence of gravity these increments remain 
superposed upon each other without alteration. Under gravity, on the other hand, each 
increment undergoes transformation in the manner characteristic of an “initial eleva- 
tion” of the same magnitude (the relevant theory is given in Lamb’s Hydrodynamics,* 
y being there measured upward). The usual Fourier integral may be omitted here, how- 
ever, since it may be verified by carrying out the integration in k that (5) can be written 





(361) 28° [cos ke dk [ F(8, te d (6 
Vv. = (| = os ka dk d e ; 
1 Nay Juno ® Jo | he 8 ) 
Now a standing wave with a surface elevation » = cos o t cos kx, where o” = gk, has 
an associated potential g(sin ot/c)e"“” cos kx (cf. Lamb, Sec. 238). Hence the potential 
for an oscillation reducing at time t = ?’ to a surface elevation »,, dé’ is, using (6), 


io = 9 sav [ o7sin o(t — t’) e~™ cos kx Qk, t”) dk, 


; (7) 
Q(k, t) = / F.(8, t) e* dg. 


The total potential ¢ is then found by integrating (7) with respect to ¢’, and adding ¢; : 
2 y : , 2 -8 ° , —ky , 
@=a+ 298] at [ o- sin o(t — #”) e-™ cos kx Q(k, ”) dk. (8) 
“0 “0 


The corresponding surface elevation is, from (3a), since 0¢,/dt = 0 at y = 0, 


OY ‘ n@o 
— “2 dt’ | cos o(t — t) cos ka Q(k, t’) dk. (9) 
Jo Jo 
It is easily verified that @ as given by (8) satisfies Laplace’s equation and the re- 
quired boundary conditions (2) and (3b). At z = 0, 0¢/dx = d¢,/dx and (2) is satisfied; 
in 0n/dt, the term obtained from the upper limit equals d¢,/dy as expressed in (6), and 
gk = a’. 

These equations represent the motion caused by an arbitrary motion of the vertical 

boundary. In the harmonic case the expressions become simpler. 
1. The harmonic case. As a special case, let the displacement of the vertical boundary 


after t = 0 be 
s = SF(y, t) = Sf(y) sin ot, (10) 
so that F,(y, t) = wf(y) cos wt. Then, replacing t/ by r = ¢t — Ut’ as the variable of inte- 


gration, (8), (9) and (7) are replaced by 


2 a a , , 
= ¢, +-agS | dr | (cos wt cos wr + sin wt sin wr) 
T /0 Jo 


(11) 


« o ‘sin or e~” cos kx K(k) dk, 


3Lamb, Hydrodynamics, Cambridge University Press, 6 ed., Sec. 238. 
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es « Ss , ; 
n = -—wS | dr | (cos wt cos wr + sin wi sin wr) 
T /0 “0 (12) 





x cos or cos kx K(k) dk, 


, (B ‘ 
K(k) = |e #(@) ag. (13) 
*o0 

These expressions represent a non-harmonic motion of the water, caused by a har- 
monic motion of the boundary after an initial state of complete rest. As t ~ © , however, 
the integrals representing the coefficients of sin wt and cos wt will approach fixed limits, 
provided f(y) satisfies the usual restrictions. Thus, the entire motion becomes more 
and more nearly harmonic. 

There must exist, therefore, an ideal solution representing exactly harmonic motion, 
in which (10) holds at all times. To obtain the corresponding amplitude functions in ¢, 
which may be written 

@ = > + vi(z, y) cos wt + ¥2(z, y) sin at, 
we carry out the integrations in 7 and take limits as t > ©. The amplitude functions thus 
found are: 
waS f° /1 — eos(o + a) 1 — cos (o — w)t\, 
“ge | (1: aS 4 (eee + 
T 


o+w ¢—-w cos wt ’ 


v(x, y) = lim 


/0 


. w9S f° sin(o +w)t , sin(o — w)t\, 
y2(z, y) lim 22° [. Set z=) dk, 


ot+w o—-w 


G =o 'e™ cos kx K(k). 


The values of the limits are obtained at once from the formulas: 
: in tz 
lim [ Ia) —— dz = zf(0) or 0 (14) 
tao va r 


according as either a < 0 < b ora and b have the same sign; 


— ff 4,1—- eek P ae 
lim | f(x) 88 dz = | f(x) “—- (15) 


too /-a J —a 


where the principal value of the last integral is taken. 
These formulas can be reduced to the familiar one, proved in many books:* 


ah 


lim | f(x) sin te dx = 0. 


t-a “a 
In (15), for example, the left-hand member can be written, because of the symmetry 
of (1 — cos tz) 


a (1 — cos tx) dz. 
H 


[ J(z) — f(—2) 


Noli 


lim 
t—@ “-a 


4Whittaker and Watson, Modern Analysis, Cambridge University Press, 3rd ed., Sec. 9.41. 
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With mild restrictions upon f(x), the first fraction here will be finite at z = 0, and the 
cos tx term will give zero in the limit; for the remainder, note that 


; dx [ dx 

2—=-— —z)—. 
[ mZ=--f 1-9 
The integral in formula (14) when a < 0 < b can be reduced in an analogous manner 


to the integral, 


| sin t2 dx/x = for [> @, 


x 


Using these formulas, the cosine terms in y, and the (o + w) term in y, are zero, 
whereas the (o — w) term in y¥, can be evaluated after changing from dk to d(a — w); 





since o = gk, dk = 20 do/g. Thus (11) becomes, using (4) and (10), 
Rae x +(B+y) 
= Iu — 7 he 
? 7 E | J(B) log x + (sp — y) - 
+ | ___ " cos kx K(k) a | COS wt (16) 
T Jo k- Ko 
+ QwS exp (—kyy) cos kyx K(ko) sin wt, 
where 
ko = w'/Q. (17) 


Nearly the same integrals occur in (12), which becomes 


2S | ss Ls a dk re). | 
— | cos kox K(k) cos wt — — sin wt | ———— cos kr K(k) dk }. (18) 
g T Jo k — ko 
The boundary conditions (2) and (3b) are again easily verified. In 0d/dy, 
—k/(k — ky) = -1 — kj/(k — ko) = —1 — w/g(k — ko); the first term serves to 
cancel the contribution from the first integral when y = 0 because of (4) and (6). 
Equations (16) and (18) represent exactly harmonic motion, but they may also be 
used at a given point as approximate expressions for the motion after a start from rest, 
provided sufficient time has elapsed. It is necessary that, as o varies past the value w 
in the integrations in (11) and (12), sin et or cos ot shall execute many periods while 
cos kx varies little from cos kyr. Consequently, a good approximation to the limits as 


t— © will be obtained. This requires that 
wt > kor = w'x/g = 2rx/dr 

where } = 2xg/w represents the wave length of small traveling waves of frequency 
w/2r. In other words, the total number of waves emitted from the boundary since the 
start, as represented by wt/2x, must greatly exceed x/\ or the number of waves included 
between the boundary and the point z. 

On the other hand, as x increases beyond the first few wave lengths from the 
boundary, the exactly harmonic motion approximates more and more closely simple 
traveling waves having the frequency of the boundary. We can show this by trans- 


forming the remaining integrals in k. 
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Putting ¢ = k + im and integrating around the first quadrant of the ¢-plane, in- 
dented at x = ky , we find that 


exp (—yf + taf) ». _ [* exp (—vk + tak) 
§ RE SD ay / pa maine 


; . [° exp (—tym — am) 
- im exp (—yk, ko) — | ne Sg as 0 

wr exp (—yky + taky) te i = i 
for y > 0 and a > 0; the principal value of the k integral is intended. The last integral 
can also be written 


“im + ky , 
- —;——-3 exp (—iym — am) dm. 


2 
Jo m + ko 
Hence, taking real parts, we obtain the formula, 


”~  k : , 
an e™* cosak dk = —r exp (—~yko) sin ak, 


m + k; 


For our present purpose, change m to k, set a = x and y = B + ying ory = 6 
in n; the factor e~™* is to be found in the integral for K, as given in (13). It is then found 
from (16), (18), (13), that 


@ = QS K(ko) exp (—koy) sin Gt — kox) 


“m cosym — kosinym _am 
+ nnd ——— @ dm. 
“0 








*+ (6+ y) 


2 bi 1c 1 x 
+p 7 oS cos wt | f(8) as| 1 log re (8 = y)° (19) 


r" k cos k(8 + y) — ko sink(B + y) - 
ee aE o-™ dk |, 


+ Je ke + ks 
Qw’ S 
9 = = | K(k) cos (wl — kox) 
tf 
(20) 
- -_ ‘° k cos Bk — ky sin Bk _y, 
— 7 sin wt i J(8) dB | - ae So e* an |. 

The first term in (19) or (20) represents a train of harmonic waves traveling out- 
ward from the moving vertical boundary. In terms of the wave length, \ = 2x/k) = 
2rq/w°, the surface amplitude of these waves is, from (20) and (13), 

4S [* 
A=" | exp (—2ny/d) fly) dy. (21) 


The second term in (19) or (20) represents a local disturbance near the wave-making 
boundary, superposed upon the waves in time quadrature with them. The resultant 
disturbance near the boundary in the exactly harmonic motion will thus have an 
amplitude in general larger than that of the distant waves. 
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At the boundary itself, the expression for the surface amplitude in the exactly har- 
monic motion becomes simpler. With x = 0, in terms of k’ = k — k, , the integral in 
(18) becomes, using (13), 


[ f(8) dB [5 —7 —— —Bk, + k’)) = - F exp (—k,)Ei(8k,) f(8) dB, 


Bi(z) = [ “a= — | oa 


Hence at x = 0, in terms of \ = 2g/w’, from (18) and (13), 
4rS 1 — ° 
ee @ cos wt + = Ei(koy) sin wt | exp (—koy) f(y) dy. (22) 


Here, as y — 0, Ei(koy) becomes logarithmically infinite; nevertheless, if f(0+-) is abso - 
lutely bounded and integrable, the integral will converge. 

The pressure on the boundary is also obtainable in closed form, for the case of 
exactly harmonic motion. It is p = p(d@/dt),-o . For the evaluation of the integral in 
(19), the following formula is easily inferred from formulas given in Bierens de Haan’s 
table of integrals: 





* k cos pk — qsin pk - , 
b= —e ™ Ki 
[SE a = - «Bion. 


where p and gq are constants. For the present case, take q = ky) , p = 8 + y. The pressure 
at depth y is thus found to be, using (13), 


p = 2w' pS cos wt exp (—koy) [ exp (—k8)f(8) dg 
Jo 


sin wt [ s(a){} to ety, > BD, 


oa) 


Qu" pi 





(23) 
— exp [—k,(8 + y)Ei[k(8 + wi} dg. 


The second term in p is a quadrature component that, on the whole, does no work on 
the boundary; the first term provides the energy carried away by the waves. 

Any motion consistent with the assumed harmonic boundary conditions can be 
resolved into the motion already described and a complementary component satisfying 
the condition of zero horizontal velocity over the vertical plane x = 0. 

2. The pure-wave case. If f(y) = exp (—k,y), the boundary moves as does a vertical 
layer of fluid particles during the passage of harmonic waves; the entire motion must 
then reduce to that of the wave train. 

Now, if f(y) = exp {—2my/d}, (21) gives at once the amplitude A = S, which 
is correct, since in deep waves the vertical and horizontal amplitudes are equal. That 
the local disturbance vanishes completely in this case is most easily seen if the corre- 
sponding term in ¢ is put into Havelock’s form. 
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The first of the following equations is easily obtained, and integration of it with 


respect to a yields the second: 


a ~az a a 
e€ (cos px — cos gx) dx = =—— —- a 
Jo a +p a +q 
r“ COS px — COS GX _a: l a’ + q 
———___ <—e dz = — log =" . 
/0 x 2 Qa of Pp 


The constant of integration in the latter formula must be zero since both members of 
the equation vanish as a >. 

Using this last formula with x = k, a = x, p = B — yandq = 8 + y to transform 
the logarithm in (19), we obtain as the coefficient of f(8)d@ in (19), 


[ k cos (8 — y)k — cos (8 + y)k ¥ k cos k(8 + y) — ko sin k(6B + | ke le 
2 k +k Pur 


Ja 
and, after consolidating, the integral in B in (19) becomes 


§ e & ° k COS ka —_ ko sin ka kz 
| f(B) dp | (k cos Bk — ky sin Bk) er +k) Y o-* dk. 


Jo Jo 
The entire expression for @ in (19) then agrees with Havelock’s expression except for 
a change of notation. Furthermore, it is easily verified, by carrying out the 6 integration, 
that if f(y) = exp (—koy), the expression just written vanishes. The right-hand member 
of (19) thus reduces in this case to the first term, representing the waves. In the same 
way the sin wt term in (20) disappears. 

Il. Case or Finire Depru 

Let it now be assumed that the liquid has only a finite depth h. Nothing essentially 
new is introduced by this limitation. The principal formulas will accordingly be written 
down with a minimum of explanation. 

In using the source method, infinite trains of line images are now required in order 
to preserve the boundary condition both on the bottom and on the free surface. The 
sources associated with an element dy of the boundary at a depth y = 8 fall into two 
series, beginning respectively at y = +8; in each series the spacing is 2h and the signs 
alternate, since reflection in the free surface reverses the sign whereas reflection in the 
bottom does not. This alternation of sign suggests use of the complex potential C log 


tan az, whose real part, if C and a are real, is (C/2) [log (cosh 2ay — cos 2ax) — log 
(cosh 2ay + cos 2ax)]. The proper periodicity is obtained if a = 2/4h; and for small 


x this expression reduces, except for an added constant, to (C/2) log (2* + y*) or to 
the potential of a line source at the origin. For one train of sources, y is replaced by 
(y — 8), for the other, by (y + 8). It is thus found that, for the motion in the absence 
of gravity, (2) and (3) representing the potential and the corresponding upward surface 
velocity are replaced by 


: FB, t) log (cosh bx — cos b(8 + y) cosh be + cos b(8 — ») dg, (24) 


, 








”” Qn J, cosh ba + cos b(6 + y) cosh bx — cos b(8 — y) ' 
ai) S *  F (8, t) sin 68 ” 

i =i = ‘osh bz "swore wet 2 
sn (2? — h cosh be i sinh” bx + sin” bg qs (25) 
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with 


The expression for v{/, is then replaced by a Fourier integral: 


fF . 
v(x, t) = [ cos ka dk v! (a, t) cos ka da. (26) 

70 “0 
The formula needed to reduce this integral can be obtained by integrating e‘** cosh 
bz(sinh* bz + sin* b8)~' dz around the upper half plane for z = x + zy. Poles occur on 





the imaginary axis where sin by = + sin b8, or y = B + nx/b, y = —B + (n + 1)x/b, 
where n 0, 1, 2, --- . The sum of the residues thus leads to a series, namely, 1 — 
exp {—knr/b} + exp {—2kr/b} --- = (1 — exp {—kn/b})~*. The real part of the 
integral is thus found to yield the formula 
r _008 kx cosh ba re == E oy 2 sinh ie | (27) 
J, sinh” bx + sin” bp 2 bsin bg om" 4 


By means of this formula, (26) with v/,(a@, t) inserted from (25) reduces to 


9) ’ h 3 iJ 
rm =" 1 coskr Qk, \dk, Qk, v= | FAB, (es + seen) dp. 
TT J “0 € 


This equation for vf, replaces (6). With the usual modifications to fit the finite depth, 
the former procedure then gives for the total velocity potential ¢’ and surface elevation 


n’ in place of (8) and (9), 


90S fr po osh k(h — 
‘=¢ 4 = | dt’ | o sin o(t — t’ cosh k(h — y) cos kx Q'(k, t’) dk, (28) 


ili ve * Jo cosh kh 
, 28 r , F , , , 
7 = — | dt | cos a(t — t’) cos kx Q’(k, t’) dk, (29) 
Tr Jo “0 
where now 
o = gk tanh kh. 
The boundary conditions are easily verified, including the new one that d¢/dy = 0 


when y = h. 

3. The harmonic case. The treatment of this case then proceeds as before. It is 
found that, in the exactly harmonic case where the displacement of the boundary is at 
all times s = Sf(y) cos wt, instead of (16) and (18), 


cosh k,(h — y) K'(k,) — 
cosh kh = tanh kh + kh sech® kh 





sin wi cos kx 





¢ =o + 2a 
(30) 





cos wt [ 
r Jo cosh kh 


cosh k(h — y) as K’'(k) dk 
— k tanh kh — k, tanh k,h 


where k, is the positive real root of the equation 


w = gk, tanh k,h (31) 
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and 
a 2 ie 2 sinh k6 
K’'(k) = ‘(8B le ss + asioh I | d ’ 32 
®= | sole" +r | 4 (32) 
, _ 2S {cos eed K’(k) 
tas g oe a we a kh + kh sech’ k,h 


(33) 





_ sin wt [" cia K'(k) dk } 
‘a Da k tanh kh — k, tanh k,h)° 

In deducing (30) from (28), an equation is first obtained that is identical with (11) 
except that ¢ and ¢, are replaced by ¢’ and ¢{ , K(k) by K’(k) as given in (32) above, 
and e ™” by cosh k(h — y)/cosh kh. With the same changes, the expressions given 
under Eq. (11) for ¥,(z, y) and y¥.(x, y) become the amplitude functions yj and y for 
¢’, such that ¢’ = of + ¥i(z, y) cos wt + ¥4(z, y) sin wt. In the further reduction, how- 
ever, k, replaces ky , and the relation o” — w” = g(k — ko) is replaced by o? — w = 
g(k tanh kh — k, tanh k,h), so that the latter expression in parentheses replaces 
(k — kp); furthermore, the relation dk = 20 da/g is replaced by dk = 20 do/g(tanh 
kh + kh sech® kh), so that at k = k, an additional factor (tanh k,h + k,h sech’ k,h) is 
introduced. The differences between (30) and (16) are thus explained. The expression 
for 7’ is then easily obtained from (3a). 

The waves at large x can be discovered as before by transforming the remaining 
integrals by contour integration. Here, the expression k tanh kh — k, tanh k,h appears in 
(30) and (33) where (k — ko) appears in (16) and (18); hence at the pole, which occurs 
here at k = k, , the following additional factor is obtained in the denominator: 


E (k tanh kh — k, tanh mt) | ; = tanh k,h + k,h sech’ k,h. 
The terms in sin wt thus combine again with the cos wt terms to represent traveling 
waves at large x. The remaining expressions are complicated, but they involve z only 
in a factor of the form exp (—k’x) where k’ is real and positive, and so represent again 
a local disturbance near the boundary. 
The surface amplitude A of the waves can be inferred from (33) with (31) and (32): 


_ 4x8 


x (1 + k,hsech k,h esch k,h)~™ 
1 


A 
(34) 





h = 2 sinh ky | 
x / f(a) exp (—ky) + exp (2k,h) + 1 7 


where k, = 27/2, in terms of the wave length A, . As h ©, this expression reverts to 


that for deep water as given in (21). 
If, on the other hand, h/), is small, so is k,y throughout the range of integration, 


and k,h sech k,h esch k,h = 1, nearly; thus, approximately, 
Y h 
A = 78 | 5G) dy. (35) 
1 0 


This is easily seen to be the correct expression for canal waves. 
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ON STEADY, LAMINAR TWO-DIMENSIONAL JETS IN COMPRESSIBLE 
VISCOUS GASES FAR BEHIND THE SLIT* 


BY 
M. Z. KRZYWOBLOCKI 
University of Illinois 


1. Introduction. Schlichting and Bickley’ solved the problem of a laminar, steady, 
two-dimensional jet in an incompressible viscous fluid flowing through a narrow slit in 
a wall and then mixing with the surrounding medium at rest. In the present paper the 
problem of such a jet in a compressible viscous fluid is solved. The following equations 
are taken into account: equations of motion, continuity, energy and state. The coeffi- 
cients of viscosity and thermal conductivity are assumed to be functions of temperature. 
In order not to obscure the problem by many items of a simple algebraic nature this 
paper presents only an outline of the method of attack. The method enables one to 
find the distribution of the velocity, density etc. far behind the slit. 

2. Basic equations. Assuming that the coefficients of viscosity 4 and heat conduc- 
tivity K are variable functions, f one obtains the following equations (equation of motion, 
continuity, state and energy): 


pel V. + grad (} V’) — (V x@)] = pF — grad p 
— #(grad »)(div V) + 2(grad u-V) V + (grad u) xo (1) 


+ 4u grad (div V) + uV’ V, 


pp + V(eV) = 0, p=RoT, (2) 

Je,p(T, + V-grad T) + pdiv V 
= J[(grad K)-(grad T) + K div (grad T)] + 4, ” 
@ = n{2V[( V-V) V] +o — 2 V-grad (div V) a 
a 


— (div V)'}, 


where @ = curl V. 








*Received Aug. 10, 1948. This paper was presented to the American Mathematical Society on 
Feb. 28, 1948. 

1H. Schlichting, Laminare Strahlausbreitung, Z. angew. Math. Mech., 13, 261-263 (1933). W. G. 
Bickley, The plane jet, The London, Edinburgh and Dublin Phil. Mag. (7) 23, 727-731 (1937). Also, 
Modern developments in fluid dynamics, ed. by S. Goldstein, Clarendon Press, Oxford, 1938. Bickley 
obtained the solution in the form of hyperbolic functions. Schlichting connected a series in ascending 
powers with an asymptotic solution, and obtained different numerical results. 

{It is quite enough to assume that both » and K are functions of the temperature only and not of 
the pressure. Both functions may be assumed to be polynomials or power series in 7 with properly 
selected coefficients (convergent series). 
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Assume steady two-dimensional flow in rectangular coordinates without external 
forces, with the z-axis in the vertical direction and with the equation of state included 
in the equation of energy. The equations may be presented in the form: 


puu, — pu,, = —p, — pu + wu, + 1/3 w,, 
(4) 
+ 4/3(uu,), — 2/3 uv, + wy, , 
puv, — 4/3 w,, = —p, — prv, + 4/3 wv, + 1/3 wu, 
(5) 
+ wu, — 2/3 wu, + (uv,), , 
up, = —plu, + v,) — p,2, (6) 
J(c,puT,,— KT,,) = —RpT(u, + v,) — Je,pT, 
(7) 
+ J((KT,), + K,T,| + 4, 
@ = wf{4/3[ur — uv, + v,] + (u, + 2,)"}, (7a) 


where subscripts denote partial differentiation. 

A further assumption is that the jet cross sections taken into consideration are 
located very far downstream behind the slit. In an application of a series expansion in 
negative powers of x this assumption enables one to neglect consistently the higher 
powers of 1/x. Consequently the expansion is an asymptotic one for large x.* 

A brief explanation of the method of solution of the set of differential equations 
will be given. The coefficients on the left-hand sides of those equations will be assumed 
so as to obtain ordinary linear differential equations of the second order. On the right- 
hand sides the results from the preceding approximations will be used. The results will 
be represented in series:** 

n “ s . 
u= dou,, v= Do y,, p= pot Dp, f=-f. + 2.7, . 
1 1 1 l 

3. Initial approximation. Assume that the gas-jet flows through a narrow slit in the 
wall of height h. The velocity distribution across the cross section of the slit being 
constant, one may write: M = hp u;,, where M denotes the rate at which the momentum 
flows across the initial section of the jet. This relation permits one to find the value wo . 
The calculations given below will present the velocity, density, and temperature pat- 
terns across the jet at some distance from the slit. The value uy will serve as an auxiliary 
parameter. 

4. First approximation. 

(a) Longitudinal velocity component. Put u, ~ x "f(n), 7 ~ yx *. Since the rate M 
must be constant in all the cross sections, one has, in the first approximation, 


*The applied method is, strictly speaking, an extension of Goldstein’s method in incompressible 
flow. See: S. Goldstein, ‘On the Two-Dimensional Steady Flow of a Viscous Fluid Behind a Solid 
Body,” Proc. Roy. Soc. of London, A, 142 (1933), p. 545-562. 

**Throughout the paper the subscript zero denotes the values referring to the conditions inside the 
container from which the gas flows out. The subscript © refers to the conditions in the undisturbed gas 


at rest outside the container. 
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M = 2p. [ uidy~ [2° *2*f*(n) da. 
Hence 2p = g. Another condition is that wu,, must be of the same degree in x as w,, , 
orp + 1 = p + 2q. These two conditions give the values: p = 1/4, g = 1/2. 

The longitudinal velocity component will always be calculated from the first equation 
of motion. In the first approximation put:* v = 0, u = wo, p = po, UU, = Up, and 
neglect the variation of the pressure, i.e., p, = 0. In this approximation only the terms 
of degree x *’* will be retained. One easily obtains 


Paotlollis ~ Hoth, = 0. (8) 


Introduce a non-dimensional coordinate 7: 


1/2 1/2 
“\4y..2 os %s y ia : ‘ 4y.2 ‘ 


(9) 
u,=A Ux ’“*f,(). 


The boundary conditions are that u, must tend to zero when x >, that u, must 
tend exponentially to zero when y —© and from the symmetry conditions that u, 
must be an even function of y, i.e., w4:, = wu, = O for y = n = 0. After some elementary 


transformation one obtains 


fl + Qnft +f, = 0. (10) 
Put 
fi = k, exp (—1) 
and obtain 


ki’ — 2nki — k, = 0. (11) 

This ordinary linear differential equation of the second order has a solution in the form 
of a series (see Appendix). The final form is 

u, = Ayua~’*k, exp (—7’), (12) 


with all the boundary conditions fulfilled.** The value of A, will be calculated from the 
condition that the rate M@ must be constant in all the cross sections. Retaining from 
k, only one term gives: 


1/2 © 1/2 
M = 4Ajus”’ pov’ (55) with [ exp (—2n’) dn = () ; (i3) 
u ite 
A, = +M' (5) a OF. (14) 
27 oP 


*See the explanation below for the selection of p.. instead of po . 
** Another solution may be obtained by putting 


wr 


fi = kis exp [—7°/2], cf, — ky. = 0 (see appendix). 


***Taking into account the second term will give the result: 


(| exp (—7') dy — 7 exp (-9) |. 


This expression will add 25% to the value of M. 
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(b) Transverse velocity component. This component will always be calculated from 
the second equation of motion. But in order to find the highest degree in x of v, one 
calculates the value of v,, from the.continuity equation for the incompressible fluid: 
U1, + Un, = 0, and obtains 


eo” 
—3/4 ; 2 ‘ 2 4 
vy = Amur * | [(1 — 49°)ky + 2nk/] exp (—7*) dy + C. (15) 

#0 
The constant of integration C must be equal to zero to fulfill the boundary conditions. 
The component v must satisfy the following conditions: it must tend to zero when 
x —o, it must tend exponentially to zero when y ~~, and from the symmetry condi- 
tions it must be an odd function of y or 7, i.e., v = 0 for y = » = O. The final result is: 


Vip = Aguox **[2nk, exp (—7’) — gual, 


: , 2 l Vo ‘as + 3/4 
gn =] kexp(—a) dn, 42 = 5 4A(22) in, 
“0 


2 Uo 


(16) 


Although not all the boundary conditions are fulfilled, it is obvious that the value found 
for v may be considered as the third approximation to v (~z~**). Hence v, = 0.* 

(c) Density. The density will always be calculated from the continuity equation. 
It must be equal to p. when t > and y >. Hence let us put p = pa + >.1 p;. On 
the left-hand side of the equation of continuity put the values u = uw, p = pr» $+ pi, 
and on the right-hand side the values u, , v; , and p = p, . Retaining the terms of degree 


z*’* one obtains: wopiz = —Potliz , OF 
pr = —pula"ts + C = —Aipor™”*k, exp(—2), C= 0." ai 


The boundary conditions for p, are identical with those for u, . 

(d) Temperature. The temperature will always be calculated from the energy equa- 
tion. In the first approximation only the terms of degree z~°’* will be retained. The 
temperature at infinity must be equal to 7. . Hence put T = T. + >>? T; . On the left- 
hand side of the energy equation put p = pa +p,,uU =~ %,7T, ~2'*,K=Ka. 
On the right-hand side put the values u, ,v,, 7 = T., K = Ka, p = po + p,. The 
boundary conditions for 7’, are similar to those for u, . One obtains 


J (C.potlol 12 — Kal iy) = —RpaoT otis . (18) 


Introduce a non-dimensional coordinate 7: 


in A Coltn\'? 
7 (ox =) fies : Beat * 
ee D cnt = CUrPo) ‘19 
a= ~ 9 a UE ~~ os) 5] ( Y) 


WaT R : 
’ ’ —1/4 —_ y ° 1/4 
Jig = Az ox h,(m), A; = (. oy A\(~ In. Je 
Je. 

*In order to verify whether the value »; = 0 fulfills the second equation of motion to the first degree 
of approximation, put into this equation on the left-hand side the values: p = pp, u = u,v =N, 
#t = feo, and on the right-hand side the values: py = 0, v = 0, u = wu, . Retain only the terms of degree 
x */*, One obtains: pouorz — (4/3)u.%yy = 0. Hence », = 0 is really a solution. 


**It was assumed that p; ~ po. 
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After substituting (19) into (18) and performing the necessary transformations one 

obtains: 
hi! + Qmhi + hy = 4[(4q? — 1k, — 2nki] exp (—71’). (20) 

Putting 

hy = Pil) exp (—7), 

one gets 

(pi’ — 2m pi — pr) exp (—7) = 4[(4n° — 1)k, — 2nkt] exp (—7’). (21) 


One complementary function of this equation is k, ;* hence put p, = kr, and obtain 


[2(k{ — ake + kir{’] exp (—7’) = Fi , (22) 
or 
1 
z [Kirt exp (—7)] = Fk, . (23) 
Hence 
- _— os 
kirt exp (—7’) = ge dn + c = Fi: ° (24) 
Jo 


As one may easily notice, after the integration is performed the expression on the right- 
hand side will consist of the following functions: the error (probability) integral and 


° 7 exp (—c7’). The constant of integration must be such as to fulfill the boundary 
condition 7; = 0 for 7 = 0. In the next step one obtains 
bt tn 
= [| EP, exp) dn +C.. (25) 
“0 


The constant of integration may be different from zero to fulfill the boundary condition 
for 7 = 0.** The function ky? may be expanded into a series by ordinary division. The 
function 7, will consist of the following functions: }°> 7°" exp (—cy’) and fi 7" erf 
(an) exp (cy’)dy. This last integral is not tabulated and must be calculated in each 
particular case by numerical methods. One term in 7, will consist of an ordinary ex- 
ponential function, which is not equal to zero for » = 0.*** Hence all the boundary 
conditions are fulfilled and the final result is 


T, A 3 YC ad exp ( _s 7); (26) 


I 


T, = A,T.x7“‘kr, exp (—B’1’). (27) 


C, will be determined from the condition that with no heat dissipation from the jet 
the total heat content (enthalpy) expressed as mass Xi = mass Xc,(7' — T'.) = 2c, 
fo uT1(p2 + pi)dy = const., in each cross section. 


*All the functions not included in the text are given in the Appendix. 
**In the present case it will be assumed that C, ¥ 0 (see below). 


***Hence 7,(0) ~ 0. But the value 7',(0) will include C,, (see below). 
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(e). Coefficients of viscosity and heat conductivity. Assuming that » can be expressed 
as pu wi+al?+b7T’+---)andKasK Koc,u, and substituting the value 7’ = 
is f (i , one obtains 


1 = to + BD, G2 , K, Koepm - (28) 


5. Second approximation. 
(a) Longitudinal velocity component. Put into the first equation of motion on the 


left-hand side the values p = pw + pi, U Ug + Uz, Us = Ure HF Use, Uyy = Utyy + 
Uoyy » & = mM, , and on the right-hand side the values p, = R(pT),, T = T. + 71, 
vy = Oand wu, ~ « '”’. In this approximation only the terms of degree 2°” will be re- 


tained: 

Up Pallas — Motlay, = BGo(T'4,), — Uopiti, — R(piT), - (29) 
All the functions have to be taken as functions of 7. The solution of this equation will 
follow precisely the procedure outlined in the calculation of 7’, . The main points are 


the following: 


Uy Ux’ f.(n), fo = k, exp [—4n'], (30 
Urpot  [fe’ + Wns + fo)] = uopox (ke! + (1 — °)k2] exp [—4n'] = F,. (31) 


One complementary function is fo, 1 + = 2,4,-- 0:7. Hence put kz = fargo : 


Unpat  (2fsig2 + fogs’) exp[—an] = F. , “ (f2193) = > Fs, 
( n l 


(32) 


v ( 


qs rif [ 2. Fs. dn + Aa], 


Because in the expression on the right-hand side of (31) there will appear terms of the 
form exp (—7’), the constant A, must be chosen so as to make g3(0) = 0, where 


g [ [ a([" Z F’,, dn + 4.) | dyn + Ag. (33) 


Since the function g, will contain terms of the form exp (—7’), which do not vanish 


( 


for 7 = 0, the constant of integration A, must be calculated and is not equal to zero.* 
The result is 

Us Ut forgo exp [—47'], (34) 
with all the boundary conditions fulfilled. 

(b) Transverse velocity component. To verify the statement made above that v, = 0, 
put into the second equation of motion on the left-hand side the values p = po + pi , 
“=z Uy + u,v Vo, M = mw, , and on the right-hand side the values p, = R(pT'), , 
u u, + U,V 0, » = mw, , and keep only the terms of degree x *”” 

SUpPaVoz — 4baVay, = O. (35) 
Hence v, = 0 is a solution. 
*The rate M with the previous value of A; and with p = p. + pi, u = UW + U2 must be constant 


in each cross section. This will give the value of Ag. 
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(c) Density. Put into the equation of continuity on the left-hand side the values 
u Uy + U2, P = Po + p; + p2, and on the right-hand side the values p = py + p, , 


u U, + Ue ,? 0. Preserving only the terms of degree x *’*, one obtains 
UoPez as — Pollo, ti Pili, ’ (36) 
p2 = —}pot *{2faige exp[—30] — Aik; exp (—2n)}, (37) 


with the constant of integration equal to zero and with all the boundary conditions 
fulfilled. 
(d) Temperature. Put into the equation of energy on the left-hand side the values 
p Bot mi +t, t= % +e, 7 = T+ 7, + T-,K = Ai; and on the right- 
hand side the values p = pn tp +ta,T=T.2.+7T, u=uH+wm,v=O4n=m. 
Retaining only the terms of degree x *” gives: 
J(C,Uppol'ss — Kaol'sy,) = JKG2(T;T;,), . a Molly 
(38) 
— JeuopiT, — Rlp.(,7; + UT.) + t,,p:T =], 


where bars over the letters denote the functions of 7. Again the procedure is identical 
with that explained above: 
T, = Tye '*h(m), he = Pz exp [—37], 
(39) 
Aon (he! + 2GqhS + he)] = Asx*"Ep!’ + (1 — 7p.) exp [-37] = G « 


With a complementary function f,,* the procedure identical to that in Sec. 5(a) gives 


fin = Bila + B,. (40) 
The constant of integration B; must be so chosen as to make 73(0) = 0, where 
rT. = | ( > B,Is; + B,) fai dj + Bs. (41) 
“oO : 1 


The constant By, is not equal to zero, hence 7,(0) # 0.** Hence: 


THY ryy a ae 9 
T, = Tox’ far. exp [— 437 ], 


(42) 
T, = Tox’ fore exp [—4}B'n'], 
with all the boundary conditions fulfilled. 
(e) Coefficients of viscosity and heat conductivity. 
Mo = Mo + BL Go(7; + T.) + G3(Ty + T>) + Gs2.T\T, + -- *}, 
(43) 


K, =  & Cy M2 - 


*The function f,, is the function f,, with 4 changed to 7. 
**B,. will be determined from the condition that 


masst = 2c, [ (po + pi + po)(u, + u2)(T, + T.) dy = const., 


) 


in each cross section. 
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6. Third approximation. 

(a) Longitudinal velocity component. Putting on the left-hand side of Eq. (4) the 
values p = po + pi + P2, U = Up + Us + Us, Us = Ure + Use + Usz, B= May Uy = 
Uryy + Usyy + Usyy , and on the right-hand side the values p, = R(pT’), , v = 0, etc., 
and letting u; ~ x *”* gives the result: 


Potlolzz ~ Hotlayy = F; : (44) 
The procedure remains unchanged: 
Us = ue “*f, , (45) 
Uppat "*(f3’ + Infi + 3f:) = Ks, 
(46) 
fs = kz exp[—3n],  uopat””*[hg’ + (2 — 4°)ka] exp [—3n'] = Ky. 
One complementary function is f;, . Put ks; = fsigs : 
9 - d vA , 
uspax™”* © (fg8) exp [31°] = fuk , 

(47) 


6 ” 6 
S393 = Lili + Cs, J3 = [ ( Y las + Chat dn + Cr. 


C, must be chosen so as to give g3(0) = 0. C; must be calculated and is not equal to 


zero. The result is: 
Uz = Ux *““*f3.93 exp [—} 7]. (48) 
(b) Transverse velocity component. Put on the left-hand side of Eq. (5) the values 
P= peta tme,’=Uwutuyt+m,v=1~2 “, uw = me, and on the right-hand 
side the values p, = R(pT), , u = wu + U2 + us, v = O. Retaining only the terms of 
degree z~’“* gives the result: 
(49) 


Spats: — 4ucVayy = Motizy ° 


One may easily show that v; = v,, , given by (16), is a solution of (49). In this case 
Ure = —Vs,* OF Urizy = —Vsyy , hence Eq. (49) changes to patovss — MoVs,, = 0. But 
v, = — J uw. dy,** or vs, = — f Usz2 dy and v3,, = — f u2,, dy. Hence after substituting 
into the simplified equation (49), i.e. po Uo f Uize dy — po f Urey, dy, differentiating with 
respect to y and integrating with respect to zx, one obtains exactly Eq. (8). But »; = 
v,, does not fulfill the boundary conditions when y ~~, viz., g,, # Ofory ~ 7-2. 


One has to find a solution of the homogeneous equation (49). To this end put 


Vs = ur “““p,, (50) 


and obtain p;’ + (3/2)np; + (9/4)p3; = 0. One distinct solution of (50) is the series 


Pa. ; another’ is 
” 
pas = Pur | pit exp [—2n"] dn + Ca . (51) 
*See Sec. 4(b). 


**The constant of integration equals 0. 
*See: E. L. Ince, Ordinary differential equations, Longmans, Green and Co., New York, 1927, p. 122. 
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An integration term by term may be applied. C; must be chosen so as to make 
(G11 + Pao) = Oas yo. 

7. Final remarks. The approximations of higher degree are only a matter of routine 
work. In the theory of incompressible fluids equations of motion and continuity are 
used. The velocity components may be derived from the stream function. But it is 
almost impossible always to satisfy all the boundary conditions for all the velocity 
components. Usually the v — component satisfies only one condition (for 7 = 0 but 
not for 7 ~~). It was shown above that in the theory of compressible flow one may 
always satisfy all the boundary conditions by use of the distinct solutions of homo- 
geneous equations. This last remark refers equally to all the dependent variables.* 


APPENDIX 


A(a). k, = 1 ae aon” ol ayn - lie + donn” ot a a2 = 1/2, a = 5/24, 
a, = 1/16, dg = 13/896, Ayo = 221/80640, etc.; 


ke = 1 + ayn + agn* ter, a, = 1/12, a, = 1/672, etc. 


4d). k=1+ DO aa’, Fu = 4{(4B°x — 1k, — 2k] exp (—B’7). 


A(e). po = BAI taT. +b +c +--:), Ga =at+ 27%. + 3To+-::, 
G.=b+3ceT.+:---, etc.; K. = Kettte . 
5(a). F. = 4[(R(oi71). + Uopitis — BGei(Uiy7's)y], b, = —1/2, b, = 1/8, 


b, = — 1/48, bs = 1/384, Dis = —7/34560, etc.; 
Fy, = Zz. F 2: ;fo. exp [—(B° + 1/2)n’], Fay = {2(As + Az) - 4{A,(1 + B) 


+ Arjn}kir,, Fue = 4(As+ Adk Kirin, Fars = (2Askirin — Arki’rs)hi , 
Foy = —Ad[(kirs), — 2B? nkyrs\(ki — 2nks), Fao = Acl2nki + (1 — 4n)kilhafar 
- exp [((—3/2)n’], A, = A,A;Rus’To(po/pa)(inch’”), 
Ay = Ai(po/pe)(inch’”), Az = A,A;Go:Topus'(inch’”). 

5(d). K = Kec, Ayo = Jelopol'n, Ge = 4{R[poltrel': + YiezTe) + thshiT' a] 
+ Jleuop:T1. — KGa(T:T)yJ— wetliy};  B, = A,AsR(Je,)~(inch™), 


B, = RT.(Jc,T,) ‘ (dimensionless), B, = A\RpoT (Jc,peT >) (inch’”), 


*This fact has been shown by the author in several papers. 
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6(a). 
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B, = A,Aspopa(inch’”), Bs; = A3G,KTopo(Kepx) (inch’”); 


In = | (4B %# — Yk, — QWeakFfn exp (—(B* + 1/2)7"] da, 


Tes = 2 | | (Bn —" 1) foro _ a(forGaal fer exp [((—1/2)(B? a 1)7'] dn, 


Tes = | (2k + (1 — 4B*9)ki)fasky exp [—(2B? — 1/2)9"] da, 


In = | ‘falki?,), — (47 — DkiF)fak, exp [—(B? + 1/2)%"] dn, 


In = — {(its)a9 aa Ankit), + (47° = Iki, + [(kirs); sis Dakss") for 
- exp [—37°/2] dn, 
Too = — | (ki — 2B-°qk,)*fax exp [—(2B™? — 1/2)3"] di; 


G; = b a 3cT’.. + tials G32 —= 2b + 6cT’.. + ahs 


FP, = BLG2:(u1,7'2 + U2,7';) + G Myla}, — Up(UizP2 + Uri) — Pallsztl2 


sii R(piT2 + poll's). ’ K; = —4P, ; far =1+ =. d;n’, 


d, = —1, d, = 1/4, d, = —1/20, d, = 1/160, diy = —1/1440, 
etc 

Iz, = —2C, | [(n° — 1) forge — n(forge)alfarki exp (—7) dn, 

Tsp = (—1/2)C, | [(4n? — Dk, — 2nki]{2fage — Aiki exp [—3n°/2]}fs 


- exp (—7) dn, 
Iss = Ax | (40? — Iki — 2nki]farfage exp (—2*) da, 


2 


2C, [ {Aiki faire exp [—(1 + B?/2)n") + (1/2) A;(2foig2 ie Atk; 


a 
~ 
I 
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- exp [—3n"/2])kir, exp [—(B* + 1 '2)n?}},f31 exp [n /2]n dn, 


en 


Ins = — | {Cale — 2nki)forre exp [—(1 + B’/2)n"] 
+ C,[( forge)» _— nf orgo\kini exp [—(B’ + I 2)n°}} S31 exp [n° /2] dn, 


In = —C | {(k — 2nk,)kir? exp [—(2B* + Da} afar exp [n°/2] dn; 


&. A,(po/px)(inch’’), C, = RpoT,(usp.) (dimensionless), 
C; = A,Go,T ups’ (inch), C, = AsGoT pu (inch), 
C; = A,A3G3;Toup= (inch**). 


6(b). psi = o(1 a > en’), e. = —9/8, e, = 63/128, 


€é, = —693/5120, ¢ -—: 6237/229376, €io = —39501/9175040, etc. 








—NOTES— 


ON TYPES OF CONVERGENCE AND ON THE BEHAVIOR OF APPROXIMATIONS 
IN THE NEIGHBORHOOD OF A MULTIPLE ROOT OF AN EQUATION* 


By E. BODEWIG (The Hague) 


1. Introduction. A chief problem of practical analysis consists in computing zeros of 
functions. For this, successive approximations are used. There are procedures in which 
the initial estimate is arbitrary and will always yield the nearby zero-point as the limit 
of a sequence of approximations. But there are only a few such procedures, the method 
of Laguerre being one. In almost all other cases the first estimate must satisfy some 
conditions in order that the sequence of approximate values will converge. How to 
obtain a useful first approximation is a special problem and has nothing to do with the 
method of approximation itself. 

Now there are two gaps in the theory of all methods of approximations. First, it is 
only proved that the method will converge; there is no possibility of comparing different 
methods of solving the same problem. Until now, one was satisfied to demonstrate the 
procedures by one or several examples and, from these special examples, general con- 
clusions about the method were drawn. Such conclusions are false, of course. 

Second, it is supposed that the convergence of the same method will be similar in 
all circumstances. The problem was not even stated as to whether the convergence 
could depend on the multiplicity of the zero which is to be approximated. On the con- 
trary, one assumed tacitly that this would not be the case, which is false again. 

Both gaps occur for the same reason, namely that until now only the notion of con- 
vergence was introduced, with no measure of the convergence. With such a measure both 
problems can be solved in principle: several approximation methods can be compared, 
and it can be decided if the character of the convergence of the same method will always 
be the same. 

2. A measure of the convergence of a sequence. A measure for the rapidity of the 
convergence of a sequence has first been given by the author.’ 

Definition. A convergent sequence 

a ae ee Sr ©. ¢ 


with the limit X is called “convergent in the degree g”’ if 


Introducing the deviations of the terms from their limit, x, — X = d, , this becomes 





lim Goes =c x0. 
a, 


To determine the degree of convergence we have, therefore, to develop d,,, in 


*Received June 21, 1948. 
1E. Bodewig, Sur la méthode de, Laguerre, Akad. Wetensch. Proc. (Amst.) 49, 911-21, especially 


912 (1946). 
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powers of d, . The lowest term cd, of this development determines the degree g and the 
coefficient c. 

The larger g is, the more rapidly the sequence will converge. Briefly we can say that, 
if k is the number of exact decimals of x, , then the number of exact decimals of z,,,.,; will 
be kg. Whatever the coefficient c may be, therefore, a method of approximation will be 
more effective, in the long run, than any other method with a lower degree of convergence. 
Only in the beginning of the sequence may the influence of the coefficient c change this 
behavior. 

For the reasons given in the introduction, the author believes the notion of degree 
of convergence to be of fundamental importance in numerical analysis. 

The reciprocal of a linear operator. We give first a general example and seek a process 
which will yield the reciprocal of a linear operator by a sequence of successive approxi- 
mations converging in degree k. 

For this purpose we consider first such a process for the reciprocal of a number N. 
Here X = 1/N, and by our definition we must have 


Gi, = Of 4- ss, 


or 
Lexi — 1/N = c(z, — 1/N)" + 

The constant c is still arbitrary. It has to be chosen in such a manner that all powers 

° , . rm: . . rk— rn 

of the unknown 1/N will be cancelled. This is done by letting ce = N*~'. Thus we have 

k E\ . . 
eae x k — . otloJem—ee tee |, where 2 am. 

2 3 

The sequence x, , x , --: will converge if 

<2; < B/N, 


and the convergence will be of degree k. 
For k = 2 the known formula 


Lar, = 2,(2 — Nz,) 


results, which converges quadratically. 

The same formula will also hold for any linear operator, since the proof did not use 
any properties other than those which are true for every linear operator. But the condi- 
tions of convergence must be established for every operator separately. 

So the formulas hold if N is a matrix. This case has been investigated by G. Schulz” 
and again by other authors independently. In a former paper’ the author has given 
another proof of the formulas slightly different from the above proof. 

3. Newton’s formula and generalizations. First we give a generalization of Newton’s 
formula. While the usual Newton formula is customarily obtained by means of Taylor’s 
theorem, we shall proceed here in another way. 





*G. Schulz, [terative Berechnung einer reziproken Matriz, Zeits. angew Math. Mech. 13, 57-59 (1933). 

3E. Bodewig, Bericht tiber die verschiedenen Methoden zur Lésung eines Systems linearer Gleichungen 
mit reellen Koeffizienten, Akad. Wetensch. Proc. (Amst.) 50, 930-941, 1104-1116, 1285-1295 (1947) and 
51, 53-64, 211-219 (1948). 
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Let f(x) have the root X of multiplicity p: 


f(x) = (x — X)”-9(z), where g(X) # 0. (1) 
Then 
f'(x) = p(x — X)”'g + («& — X)’o’ 
or 
pers heb 2) 
Near x = X the first term on the right is very large, while the second one is “finite’’. 


So with very good approximation, 


f(x) _ sp 
f(x) 2—X’ 
that is 
’ f(x) 
X =2r- p> 
, Py (x) 


Thus the general Newton’s formula in the neighborhood of a root of multiplicity p is 


f (x2) 
La+1 = 2, — P f"(z,) ° (3) 
More generally, we set 


(rs) 3" 


Fn , = i a 7 
, Ft (an) 
and determine the character of the convergence of the sequence 
ry 9 ro 9 3 be ee A i tie —> X, 


which we suppose converges, the conditions on x, being unknown. Developing g(x) at 


¢ = X we obtain from (1) that 


f (an) = (at, — XY’g(X) + (2, — X)?""'9(X) + - > 


= d?g(X) + de*'g'(X) + +: 


or by (3’), 
a a g'(X) » 
la @ 2 ~- -4*6¢ 35S Gt 
Pp p g(X) 
Subtracting X from both sides, we get 
"(X) » 
dusy = (1 - )a,+ $M ay... 
p p g(X) 
Here the linear term of d, will be cancelled only if a = p. So we have the following 


results. 
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Result 1. In the neighborhood of a root of multiplicity p formula (3’) will yield a 
sequence which converges quadratically to the root (if it converges at all) only for a = p. 
The convergence is better the larger p is. 

2. For a + p, however, the sequence is only linearly convergent; the smaller the 
quotient a/p is, the worse the convergence. So Newton’s formula will not give better 
convergence in this case than, for instance, the rule of false position, which is often 
simpler. 

3. The usual Newton formula (having a = 1) will therefore converge quadratically 
only near a simple root, but linearly near a multiple root. 

4, If it is known that the root is multiple (for instance by observing that f’(x,) is 
small) but the exact multiplicity is not known, we put a = 2 or perhaps a = 3. 

A formula which always yields quadratic convergence, without requiring knowledge of 
the multiplicity of the roots, is obtained by choosing a function which has all the roots 
of f(x), with each of them simple. This is f/f’. Then (3’) with a = 1 will yield 


Ln+1 = Ly — | (4) 
n+1 n f 2 < ff 44 
where on the right side the argument z, has been omitted. The calculation is simpler 
if we put (4) in the form 


Ins =%e— pe with A=ff"/f’. (4’) 
More generally we put 


Inti = IX — wil, (5) 


and determine the character of the convergence for the case that 
f(x) = (@@ — X)’g(z). 

Proceeding as above, we obtain 

daai1 = C14, + Cod, + °°, 
where 
c= —(p—1)(a—1)/(p+a-—ap), c= —(a+ap— p)g'(X)/p(p + a — ap)’g(X). 
So c, vanishes in the cases a = 1 and p = 1. The first case yields formula (4). In the 
second case also c, will vanish if a = 1/2. So we have the results below. 

Result 1. Formula (4) yields a convergence which is always quadratic. Since the 
factor c, is equal to g’(X)/pg(X), the convergence is least good in the case of a simple 
root; the greater the multiplicity p of the root, however, the better the convergence. 

2. The quadratic convergence of (4) for a root of multiplicity p is worse than that 
of (3) if p ¥ 1 and is the same if p = 1. 

3. The simple formula (3) is therefore always preferable to (4), though only applicable 
if p is known. 

4. In the case of a simple root, formula (5), with a arbitrary, will always yield quad- 
ratic convergence. It is, therefore, no better than the usual Newton formula. The 

«Newton formula, moreover, is preferable as it does not contain f”. 
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5. Formula (5) with a = 1/2 yields cubic convergence in the case of a simple root, 
with 


,- -12@ 4 [gaol 
3 2 g(X) g(X) ; 


but only linear convergence with c, = (p — 1)/(p + 1) in the case of a multiple root. 
Hence it is no better a formula than the simple Newton formula would be. 

6. A convergence which is always cubical, is obtained by taking the average of the 
results of formulas (3) and (4). For the coefficients c, are absolutely equal, but of opposite 
sign in both formulae. 

The simplified Newton method. It is known that the usual Newton formula z,,, = 
Ln, — f(x,)/f'(x,) is often simplified to 





f(z,) 
f(x) 


with constant denominator. Some authors (for instance Whittaker-Robinson) consider 
formula (6) to be nearly as effective as the usual one. This is, however, not the case. 
For it is easily seen that every sequence 


(6) 


Ln+1 = Ln — 


2 


Inti = In — Cf (Ln) 


yields d,4, = d, — q(x.) = d, — cdrg(X) + --- , and so always converges linearly 
for every p. i 

The same result Has been given by Ostrowski.* The present proof is much simpler. 
It must be added, however, that for Ostrowski the chief point is to give the conditions 
for convergence. 

Example. To illustrate the difference between quadratic and cubic convergence we 
take an example of H. S. Wall.*° The square root of b is calculated by means of Newton’s 
formula and by means of (5). Taking f(z) = 2” — b we obtain the two formulas: 


1/20”, 


Insi = 3(2, + B/z,) with Co 


Insic= (x2 + 3b)x,/(8227 +b) with cs = 1/4b. 


Putting b = 2 and x, = 1 we obtain the sequences where all figures are correct and only 
the last ones are slightly different: 


quadratic cubic 
convergence convergence 
2, 1 1 
Le 1.5 1.40 
Z3 1.417 1.4142132 
Xs 1.414216 1.414213562373095048796 


Ls 1.414213562375 








‘A. Ostrowski, Uber eine Modifikation des Newtonschen Naherungsverfahrens, Akad. Nauk. SSSR 
(Gruz, Fil Trudy) 2, 241-249 (1937). 
5H. S. Wall, A modification of Newton’s method, Amer. Math. Monthly 55, 90-91 (1948). 
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Starting from x, = 10 we have the sequences: 
quadratic cubie 
convergence convergence 
be 10 10 
Le 5 30 
Ls 2.7 1.6 
Xs ba 1.415 
Ts 1.44 1 .414213564 


4, Laguerre’s method for algebraic equations and related formulas. While Newton’s 
method and the generalized methods above are applicable to all equations and hold 
for real and complex roots, Laguerre” has given an elegant method which applies only 
to algebraic equations having ali real roots. The sequence is defined by 


nf (2) 


~ f(a.) + [H(2x,)]'” 


Levi = X 


where n is the degree of the equation f(z) = 0 and H(x) = (n — 1)?f” — n(n — I)ff” = 
Hesse’s covariant of f. The sign of the square root can be taken as either + or —, for 
the method approximates two roots simultaneously. 

The method has been analyzed in an earlier paper of the author." The features of 
the method are: 

1. Laguerre proved that the method always converges, i.e. 2, is arbitrary (in con- 


trast to Newton’s method, for example, wherein x, must satisfy certain conditions). 
2. Laguerre proved that each sign of the square root yields a convergent sequence: 


Ee, v r ry 
Deg Be 5 Me 5 88 5 Be ttt OS 


‘Ad ‘Ad ‘7 ree 
Denaihe oe 9° ** He p32 Sa 


where X’ is the nearest root lying to the left of x, and X” the nearest root lying to the 
right of x, . This will always hold if the straight line is considered closed at infinity. 

3. The author has proved that the convergence is cubic in the case of a simple root, 
but only linear in the case of a multiple root. 

4. In order to obtain cubic convergence in the case of a root of multiplicity p, the 
author proved that in the above formula, H(x,) is to be replaced by aH(x,) where 
a= (n — p)/p(n — 1). 

5. Later van der Corput’ proved that if the equation has at least one root of multi- 
plicity 2p, the formula mentioned which replaces H(x) by aH(x) will approximate 
only the two roots of multiplicity =p lying to the left and to the right of x, , thus skip- 
ping the roots of multiplicity <p. In the case of a root of multiplicity p the convergence 
is cubic, whereas in the case of a root of multiplicity >p it is only linear. 

5. The general convergence of degree n. In the preceding some important examples 


6. N. Laguerre, Sur une méthode pour obtenir par approximation les racines d’une équation algébrique 
qui a toutes ses racines réelles, Oeuvres de Laguerre 1, 87-103 (1898). 

7See footnote 1. 

J. G. van der Corput, Sur l’approximation de Laguerre des racines d’une équation algébrique qui a 
toutes ses racines réelles, Akad. Wetensch. Proc. (Amst.) 49, 922-929 (1946). 
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have been given showing the difference in convergence near a simple and a multiple 


root. This will now be generalized. 
For this purpose we suppose we have a process F(x) yielding to every value 2, a 
(7) 


better value 2,4, : 

Saex = F(t). 

We are not concerned here with the conditions for convergence itself of the sequence, 

but merely with the character of the convergence. Now we seek the general function 
F(x) which yields a sequence converging in degree n. Then 

Yen. — X = C(x, — X)* + --- (8) 

d‘/dx’ is the ith derivate, 

--(n— t+ Il@ — 2)” +: 


Now when D’ 
= D's,,., = cen(n — 1) 


D'F (x,) 
Putting 2, = X, the right side will vanish if 7 < n. So 
D‘F(X) = 0 for ¢=1,2---,2-1. (9) 
This is the modified condition that the sequence 2, , 22 , converges in degree n. 
3efore determining F(x) we conclude from (8) that 
X = lim z,, = lim F(z,) = F(lim 2,) = F(X), 
that is 
X = F(X). (10) 
This condition holds for every root X of f(x) = 0. Therefore 
F(x) = x + f(x)g(x), (11) 
where g(x) is arbitrary. 
Now if (9) is to hold for 7 = 1, the sequence x, must converge quadratically. There- 
fore, since 
0 = DF(X) = 14 f'(X)g(X), 


g(X) = -1/f(X) or — g(x) = —1/f'(x) + f(z) giz) 
11) becomes 
F(a) = x — f(x)/f'(a) + fw) g(a). (12) 


[It will yield a quadratically convergent sequence for an arbitrary g,(2). 


and the function ( 
Now if (9) has to hold for 7 = 2, that is if the sequence has to converge cubically, 


0 = D?F(X) = —D*(f/f’) + D?(f?9:). 
Here f’g, is differentiated by Leibnitz’ rule, and in the result x is replaced by X. Then 


all terms will vanish, except the last one: g,D?(f?) = g,(X)-2f(X). On the other hand 


we obtain 


= X, if for brevity we set r = 1/f’, 


since f(X) = 0, we have for x 
= 2Df-Dr + rD*f = —2f'f"r + f'r = —f"'r. 


D?(f/f’) = D*(fr) = 
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Therefore 
0 = D'F(X) = f"(X)-r(X) + 2f"(X)g.(X), 
or 
g(x) = —1/2 f"(a)r°(x) + f(x) go(z) 
and 


SO ae en ¢ °2 e773 3) (12) 
F(z) =2z-—fr V/2f f'r + fo, (13) 
where g2(z) is arbitrary. 
To obtain generally a sequence converging in degree n we generalize (12) and (13) by 
means of the operator 
P = r-D, where r = 1/f'(2) and D = d/dx 
and write 
F(z) = F,(x) — f*(x)g,(z), (14) 


where 


n—-1 
( 1) go~! pet, (14a) 


+ / l p2 I pS 2 
Fig) =~ fetal Pr- ef -Prt+ --- +o hy! 


We have to prove that this is the expression we obtain if we proceed from (13) to 
sequences which converge in degree 4, 5, --- , n. Since the form of F, and F; is identical 
with (12) and (13), respectively, we proceed by induction and assume that F,(z) is 
the true expression; we have to prove that then F,,,(x) is the expression yielding a 
sequence which converges in degree n + 1. 

Now since 


D(fr) =r Df +f Dr=1-+ ff’ -Pr, 


DF(z) = 1 — (1 + ff’-Pr) + (ff’-Pr + 4f?f’-P’r) — 48 f'*f’-P’r 


“3 r (—1 —2 n—2 a-le a-l 
+ f*f'-Py) + --- +e — YO Poe + ff Ps 


(n — 1)! 


= a ig f’ F*'*. 
Then since 
D"(f(X)]" = mIlf’(X)]", 
it’follows that 
D'F,(X) = D™"DF,(X) = (—1)" Uf" (Xf (X) «Pr. 
Furthermore, 


D"{(f(X)T"gn(X)} = nlf’ (X)on(X), 
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that is 
r ] n—1 
GAX) = a ye q.e.d. 


Although (14) and (14a) have been constructed to yield a sequence converging in 
degree 7 (if it converges at all), we have to investigate what happens when the root X 
is multiple. In this case the equation 0 = DF(X) standing between Eqs. (11) and (12) 
can no longer be satisfied by any function g(x) whatever. For we must have 


0 = DF(X) = 1+ f’(X)g(X) = 1, 


since f’(X) = O and X is multiple. So not even a quadratic convergence can be obtained, 
that is, the convergence is linear. 

The formula (14) can be given another form which originally dates from Euler.° 
He gave the following generalization of the method of Newton about which the author 
has reported in a former paper.’® 

For the function y = f(x) a value X has to be determined such that y = 0. Using 
the inverse function z = u(y), X is given by X = u(0). But wu is unknown, and of X 
only an approximate value a is known. Thus, f(a) is known. Now, for A = f(a) Euler 
develops u(0) in a Taylor series near A: 


X = u(A — A) = u(A) — Au'(A) + 1/2A*u'""(A) — 1/6 A®u’"(A) +--+.) (15) 
Here the derivates of u for the argument A are known, for 
u(A) = a, u'(A) = dx/dy = 1/f’(a) = r(a), 
u"(A) = —f’’(a)r*(a), u'’’"(A) = [3f’"(a) — f’(a)f’”’(a)]r°(a), --- 


Setting 


u™ (y) = (—1)""'U,,(z)r”" (2), 
we have the recurrence formula for U, : 
Ones = Gan — 12)f"U, — f°. (16) 


Now the first terms of (15) are identical with the corresponding terms of (14a). 
This must be proved true for all the terms. Thus we proceed again by induction and 
suppose that (15) holds for n and prove that it will hold for n + 1. Let therefore 


Pr = ua’ (y), then 
P*r = PP’r = Pu™ =r-Du™ = (—1)""'7[ Ur" — (Qn — 1)f"’U "| 


= (—1)°U,,,7""* = u"*(y), q.e.d. 


So we have the results following. 

Result 1. The developments (14) and (15) are identical. 

2. The convergence of the sequence 2, , 22 , --+ yielded by the function F(z) of 
(14) and (14a) is of degree n if the root X to which it converges is simple, but only 
linear if the root is multiple. 

L. Euler, Institutiones Calculi Differentialis, II, Cap. IX.—Opera Omnia, ser. I, vol. X, p. 422-455. 

0). Bodewig, Uber das Eulersche Verfahren zur Auflésung numerischer Gleichungen, Commentarii 
Math. Helvetici 8, 1-4 (1935). 
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ON THE COMPRESSION OF A SHORT CYLINDER 
BETWEEN ROUGH END-BLOCKS* 
By F. EDELMAN (Brown University) 


1. In this paper we consider non-uniform elastic compression between parallel end- 
blocks of a cylindrical test specimen with circular cross-section. It is assumed that there 
is sufficient friction between the end-blocks and the end faces of the specimen to prevent 
slippage. Thus, points on the end faces of the cylinder undergo only axial displacements. 

It is the purpose of this paper to determine the relation between the force F acting 
on the specimen, the compression 2a, and Young’s modulus £ of the material, assuming 
Poisson’s ratio »y to be known. If lateral expansion were not obstructed by the friction 
at the end surfaces, these quantities would be connected by the relation EH = Fh/Aa, 
where 2h denotes the height of the specimen, b its radius and A its cross-sectional area. 
If there is friction at the end surfaces, a correction must be applied to this simple formula: 
the derivation of this correction will now be sketched. 

The above problem is solved to a good approximation by finding upper and lower 
bounds for the strain energy of the specimen, using the method of Prager and Synge.’ 

It is advantageous to resolve the original problem, which we shall denote by Pp , 
into two simpler problems P, and P, where P, is the problem of simple compression 
and P is that problem which will yield the initial problem P, when superimposed on P, . 

In the natural cylindrical coordinates, the original problem P, involves the boundary 


conditions 


u,(r, +h) = 0, u(r, +h) = ¥a, 
(1) 


I 
= 


T16:2) =, T ,(b, 2) 


where u, and u, are the displacements in the r and z directions respectively and 7’, , 
T, are the normal and tangential components of stress transmitted through the surface. 

Thus the problem P, with which we shall deal from now on, has the boundary 
conditions 


u(r, +h) = —var, u(r, +h) = 0, 
T,.(b,2) =0, T(b,2) =0. 


The analysis used in the following is quite similar to that given in a previous paper’ 
in which the case of plane strain was considered, except that here we have cylindrical 
coordinates and different boundary conditions. The similarity, however, permits the 
procedure and results alone to be indicated here. 

2. Since the bar possesses rotational symmetry, the displacements u, and u, are 
functions of r and z only and not of @, and there is no circumferential displacement, 


i.e. Ue = 0u/00 = 0. 
*Received Sept. 22, 1948. 
1W. Prager and J. L. Synge, Approximations in elasticity based on the concept of function space, 


Q. Appl. Math. 5, 241-269 (1947). 
2H. J. Greenberg and R. Truell, On a problem in plane strain, Q. Appl. Math. 6, 53 (1948). 
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Hooke’s law then has the form 


i 
Cre = 7G [En — (Ew + E.,)), | 
] 


| 
.=@§= E E,. — v(E,, + Ee)], 
, 

















ee : , 3) 
Coo = E [Koo - W(E,, -+- E..)}, ( 
l 
é, = Ke ok E, | 
E 
Co = C9 = O 
and 
E . (1 er, + v(eos + e.:)] 
Sp St eee — ve,, v(e, €.:)|, 
(1 + »)(1 — 2p) Pon e 
| 
E 2 (1 e., + v(e,, + eos)] | 
ee ere ee — Wess v(e,, €o0)|5 
(1 + »)(1 — 2) . v1 
(4) 
; 
Eu = 7 ———~ [(1 — v)e v(e,, 2], 
(1 + »)(1 — 2v) | ee ee 
, E , : | 
Log = i+.” . En = 4 = ©. | 
The strains are defined in terms of the displacements by 
\ 
Ou, Ou, | 
C,, = ; e,, = —, 
or 0z 
(5) 
u, l (Ou, aus) 
The equilibrium equations reduce to 
9 : “a ? F a : 9 : . 
: (rE,,) + : (rE,.) = Eo, © (rE,.) + — (rE,,) = 0. (6) 
or 0z or Oz 
The Eqs. (6) are satisfied automatically if we choose our stresses to be 
‘ . l dy 
E,, = ¥, E,, = —->, 
r r or 
(7) 
’ oy rie) 1 dg 
En = = —s, ..2=->, 
” or T Oz” 7 7 O82 
where w and ¢ are two arbitrary functions, called the stress functions. 


Following the Prager-Synge notation and development we have the following 


boundary conditions: 
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a) On §*, the completely associated state, 

u*(r, +h) = —var, ut(r, +h) = 0; (8) 
b) On S{ , the homogeneous associated states, 

pu(r, +h) = 0, pui(r, +h) = 0; (9) 
c) On 8%’, the complimentary states, 

f7i(b, z) = 0, g74(b, 2) = 0, (10) 


where the stresses ,/{/ must satisfy the equilibrium conditions (6). 

3. If S,, S. and S are the “natural” states of the problems P, , P, and P respectively, 
then it can be seen at once that the strain energies can be superimposed by simple 
addition. 

Thus 


Sass (11) 


since the boundary conditions of P, and P are such that the scalar product S,-S = 
J TSu, ds vanishes. Finding bounds on S> therefore reduces to finding bounds on S’ 
since S? is known to be 2ra°L. 

Prager-Synge give the following inequality for S’: 


> (s*- 2 < 8 < S*” — > (S*- I)” (12) 


q=l p=l 


; , 


where If , I,’ are obtained by orthonormalizing the sequences S$} , S{’ respectively. As 


in footnote 2, this reduces to 
9 PF |) F 
—_—— Ses -— so 13 
r(U, + 2)a — — r(L, + 2) a’ (13) 
where 
I ° l . 
i, = — 35 (S*-77’)’, U, = —z| S* —- S*-1/)’ |. (14) 
ra E 2 a) "  mwak p> ( ‘ 

4. Suitable displacement functions, satisfying appropriate boundary conditions and 
obvious symmetry conditions, are chosen to obtain the states S* and S/ . Stress functions 
gy and y are selected for the computation of the states $/’. The orthonormalization is 
carried out using the Gram-Schmidt procedure. 

The bounds U,, and L, are computed at each iteration and after five steps we obtain 
the inequality 


F ? _ F 
0.29958 — < E < 0.30371 —, 
a a 
where we have now set b = h = 1 andv = 1/3. 
Averaging, we have 
ae ‘ 
E = 0.3017 me (15) 


In terms of the mean stress o = F/xb’ and the mean strain e = a/h, we obtain 


E = 0.9477 ©. (16) 
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Let us now perform the test in which the specimen is in a state of compression, such 


as in problem P, 


5 Fh 
S= Aa 


og 


€ 


, and measure the applied force F and the compression 2a. The formula 


(17) 


no longer yields the true value of EF since this formula applies exactly only to the test 
in which there is no friction at the end-blocks, i.e. simple compression. However the 
true value of Z is given by Eq. (16) and the relation between E and E is found by com- 
parison of (16) and (17), namely 

0.9477 E. 


E= (18) 


This gives us the correction factor to be applied to Eq. (17) to yield the true value 
of Young’s modulus Z. The numerical coefficient of E in (18) is in error by less than 
0.7%. 

The displacements and strains of the completely associated state and the homo- 
geneous associated states are given in Table I below. The stress functions and stresses 
of the complementary states are shown in Table II. 


TABLE I. The Associated States 
















































































- “ 7 | . a e, | Che Cis Cre 
s; y -“ 7 a —2)z 0 | 0 —(1 — 32’) 0 
| Sf wa 2) | 0 2(1 — 2°) | 2°(1 — 2’) 0 ra(1 — 22”) 
_ 0 | 70 — 2°)z 0 | 0 —r'(1 — 32’) | —re(1 — 2’) 
| S: R A — | 0 Qrz’(1 — 2°) | re"(1 — 2’) 0 r'2(1 — 22”) 
| s/ | 0 | ra —2)z 0 0 —r(1 — 32”) | —32(1 — 2) 
E j pare | 9 —4z' —42* 0 —3rz' 
oO TABLE II. The Complementary States. 
| | g | v E;, Et Eu EN 
Sy | hz'r*(1 — 1°) | 0 0 62r°(1 — r*) | 24(2r? — 1) | 22*r(1 — 7°) 
! "2 | oo | zr(l — 7’) | 2(1 —7’) | 21 — 37’) 0 0 
ps | ¥ 0 0 0 - 0 
| a | 0 | 2r'(1 — 7°) | 2r(1 — 2°) | 2e’r(1 — 2r’) a 0 0 
| sy | # | 0 0 | 0 =r 0 
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RIGOROUS SOLUTION OF A DIFFERENTIAL EQUATION IN 
SOIL MECHANICS* 


By R. GRAN OLSSON (Norwegian Institute of Technology, Trondheim) 


In one of his earlier books K. Terzaghi’ published a differential equation describing 
approximately the progress of consolidation in a sediment which is being deposited at 
a constant rate q per unit of time. In connection with his solution of this equation 
Terzaghi pointed out that it does not satisfy the boundary conditions for time ¢ = 0. 
In the following paragraphs the writer presents a solution which satisfies the boundary 
conditions for any time, t = 0 included. 

Let a,[gm~'cm*] be-the coefficient of compressibility, C a constant of integration, 
c [sec] a time constant, k[em sec”'] the coefficient of permeability (Darcy’s coefficient), 
qigm cm~’ sec] the quantity of sedimentation per unit of area and time, ¢ [sec] the time, 
y:lgm cm~*] the unit weight, y,,[gm cm™*] the unit weight of water = 1 [gm cm™], 
y [gm cm~*] the submerged unit weight (=y: — yw), ¢ the time function, z = (c/t)'” a 
dimensionless independent variable, and ¢(z) the error integral. 

Using these symbols, which are identical with those used by Terzaghi, except that 
Y. represents the unit weight of water, the differential equation” assumes the form 


dq , 3 (! “) c ( 
4. < aS. J 1 

at’ #\27 3/8 = # 

where c = 3yk/y, ag’ is a constant with the dimension time. 
With 
z= (c/t)' 
then t= cz° 
lg 2 dt 

and lt = —2cz” dz, > = — — 
= > ce dt 2c dz ’ 


substituting these values into Eq. (1), we obtain the equation 


de (3 = o. , 
li (: + 22)) — 22. (la) 


This equation is solved in the usual way by first finding the complementary function 


from the homogenous equation 


Integrating we obtain 


¢(z) = Cz exp (z°) 


where C is an arbitrary constant. 
*Received Oct. 27, 1948. 
1K, Terzaghi, Erdbaumechanik auf bodenphysikalischer Grundlage, Franz Deuticke, Leipzig and 


Vienna, 1924, p. 175. 
2K. Terzaghi, op. cit., p. 175, Eq. (119). 
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The general solution is given as 


¢(z) = C(z)z” exp (2°), 
where C(z) is an arbitrary function to be determined by (la). Differentiating we get 
dE _ eeyg® L (202 4 eV) 2 
= [C’(z)z” + (32° + 22°)C(z)] exp (2), 
az 
which inserted into (1a) yields, 
C’(z) = —22°° exp (—2’). 
By integration we finally get 
C(z) = -—2 | exp (—z)z ~~ dz. 
To find the integral we integrate by parts: 
F — exp (—2z°) 1 
| exp (—z)z ° dz = — =P - —2 | exp (—z°) dz, 


and thus 


C(z) = 2 exp (—2) +4 | exp (—2°) dz + C, , 


~ 


where C, is again an arbitrary constant, which is determined by the boundary condition 


fort = 0. 
a ° ° . 3 
It is convenient to introduce the error integral as a known (and tabulated) function’, 
j \ 2 l i, 
exp (—2°) dz = 5 3 (2). 
Thus we get for C(z) the expression 
9 
Y & 2 ¢ /2 y 
C(z) = = exp (—2’) + 2r'’o(z) + C, 


and for the complete solution of the Eq. (1), 


if 
- 


“(z) = 22° exp | exp (—2°) + 2'“o(z) + = c, | 


The constant C, is determined from the boundary condition that for t = 0, ¢ = {> = 1. 
When t > 0, z >=, we resort to the asymptotic development of ¢(z) in a semiconvergent 


series 
exp (—z’) ( 1 1-3 1-3-5 ) 
z)=1|1- —_ {-—-=-=s3 ———- a nes — 9-8 
22° + (2z°)° (22”)° T 
and get for ¢(z), 


' re yo > 1 > l 
‘(z) = 22° exp ©: exp (—2") + #'” — — exp (-2(1 ~ 2a 





1-3 = 


E. Jahnke and F. Emde, Tables of functions, Dover Publications, New York, 1945, p. 24. 
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By putting C, = —2z2'” we get ¢ = 1 when z >~, since both constant terms and z™* 
exp (z”) and —z”* exp (z’) cancel. We finally obtain 

oid oF a pniexp(—¢) 1:3 . 3 

(2) = 2z exp (¢)( SBE) ~~ exp (—z) t+ —---J=1- a7 se 


This satisfies exactly the boundary condition ¢ = 1 fort = 0. The diagram shown in 
Fig. 1 gives ¢ as a function of the non-dimensional “time” t/c. It is to be observed that 





10- 


08- 

















0 5 ly 10 


for small values of t/e we must use very accurately tabulated values of ¢(z) because 
the function ¢ appears as the difference between two nearly equal quantities.* In 
general, the solution of the differential equation satisfying exactly the boundary con- 
dition ¢ = 1 for ¢ = 0 will have the form 


¢(z) = 22° exp (z*)[z"’ exp (—2”) + 2'o(z) — 2°”). 


t 
4J. Burgess, On the definite integral (2/x'’”) [ exp (—?f) dt with extended tables., Trans. Roy. 
/0 


Soc. Edingurgh 39 (Part II), (1898) 
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This solution may be given a different form by introducing the derivative of the error 
function 


P 2 lie eas 
$'(2) = —7 exp (—z) or exp(?) = a7 [¢’(z)]. 


Thus we get 
f(z) = 4[g’(2)]"2"[o(@) + 227'¢'(z) — 1) 


as the most elegant form in which the solution may be written. 
Besides ¢(z) the derivative d¢/dt is also of interest. From the Eq. (la) we get 


“ = 42°[’(z)]~"[o(z) + 227'o’(2) — 1](B27' + 22) — 2z. 


Further we have 


d 


“>; 


dg dgdz _ = 
dt dzdt 2c 


Qa 


ae 
and finally 


ihe 2. 6 - _ 7 : 
z ites = [6’(z)]~"[o(2) + 227*6’(z) — 1](B2* + 22) +e | 


The diagram in Fig. 1 gives d¢/dt as a function of the independent variable t/c. 


APPENDIX. When the sedimentation process is finished, the time factor increases 
according to the equation [1, p. 176, Eq. (122)] 


(1 — §) exp (2kt/aywhi) = (1 — $1) exp (2kt,/ayhi) (1) 


or by introducing the ‘time constant” 


__ 3k 
c  y~aq @) 
into the Eq. (1) we get 
(1 — §) exp (2ct/3t;) = (1 — &,) exp (2c/3¢,). (1b) 


The time is so determined, that ¢ is equal to zero at the beginning of the sedimentation 
process and ¢ = ¢, at its end. We obtain the value of ¢, from the earlier solution of the 
differential equation, i.e., 


i= so'(e"2\| ed + 53 9) - |, (3) 


by introducing z, = (c’/t,)’”” into the solution (3). 

Equation (1b) determines a system of curves starting from a point (c/t, , £,) on 
the curve ¢ and having the horizontal line ¢ = 1 as an asymptote. Thus we get the 
system of curves represented in Fig. 2. It can easily be proved that there exists only 
one set of such curves according to the Eqs. (1) and (1b), so that the curves in Fig. 2 








342 NOTES [Vol. VII, No. 3 
are a general representation of the consolidation problem after finishing of the sedi- 


mentation proce In Fig. 2 the upper curves correspond to the greater values of k/a 
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as valid for sand and other permeable materials, whereas the lower curves correspond 
to clay and similar earth masses with smaller values of the coefficient of permeability. 


NOTE ON THE PROBLEM OF TWISTING OF A CIRCULAR RING SECTOR* 
By ERIC REISSNER (Massachusetts Institute of Technology) 


The problem of twisting of a circular ring sector is of some interest in connection 
with the calculation of stresses and deformations in close-coiled helical springs. To be 
considered is a ring-sector under the action of two equal and opposite forces P along 
the axis through the center of the ring and perpendicular to the plane of the ring (Fig. 1). 
A formulation of the problem and an outline of results by O. Géhner for sectors of solid 
circular and rectangular cross section may be found on pp. 355-361 in Theory of elasticity 
by 8. Timoshenko. 

The purpose of the present note is to obtain explicit results for the twisting of ring 
sectors of hollow cross sections, with thin walls. Formulas will be obtained which have 
the same meaning for the present problem as R. Bredt’s formulas have for the problem 
of St. Venant torsion of cylindrical rods. 

The problem may be considered as one of the membrane theory of thin shells of 


*Received Dec. 23, 1948. Work on this note was supported by the Office of Naval Research under 


Contract N5ori-07834. 
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PR 

















revolution, with multi-valued expressions for displacements. A direct solution, without 
recourse to general shell theory may be obtained as follows. The assumptions of the 
twisting theory of ring sectors are equivalent to requiring that all stress resultants of 
the membrane theory vanish with the exception of the shear stress resultant S acting 
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over the cross section of the ring sector. (Fig. 2). As all stresses possess rotational sym- 
metry the resultant S satisfies the following equilibrium equation 
dr _ 


- 
a?) +. g 


ds ds 


0, (1) 
where s is the are length measured along the center line of the tube wall and r is the 
distance of the points of this line from the axis of the ring sector. From (1) it follows that 

rs = c. (2) 
Note that the largest value of the shear stress resultant S occurs at the point nearest 


the axis. From the relation 


P= £ Ssin¢dds (3) 


the value of the constant C is obtained: 
1 
C=P/$S. (4) 


Equation (2) may thus be written in the following form 
P 
FUP a 

Equation (5) is the counterpart of Bredt’s formula S = 7'/2A of torsion theory for 
the shear stress resultant S in terms of the twisting couple 7’. 


For the determination of the deformation of the ring sector we use the stress-strain 
relation 


S= (5) 


S = Gty (6) 


where G is the modulus of rigidity, ¢ the wall thickness and y the appropriate shearing 
strain component. In terms of the components of strain in cylindridal coordinates 
Yo, and ys, we have 

Y = Yor COS + Yo, Sin >. (7) 


In view of the fact that S is the only non-vanishing stress resultant, all other com- 
ponents of strain for cylindrical coordinates vanish, 


é=@=—€& =, = 0. (8) 

Components of displacement in the radial, circumferential and axial directions which 
are compatible with (8) are* 

U = 0, V = V(r, 2), W = ko. (9) 


From (9) the remaining components of strain are obtained in the form 


=r2(Y)+1U_,2(27) 
Yor =" or \r r 00.360 @r\r/’ 


aV , 1W_aV kk 


dr * r 00 + @z r- 


(10) 





*These expressions hold without any assumption concerning the form of the cross section of the ring 
sector. 
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Combination of (10) and (7) gives 


0(V\dr OV dz k dz 
ati or (") ds az ds r ds (1 1) 
or 
y_d (¥) k dz 
rT - ds Tr . r° ds (12) 


In Eq. (12) one takes y in terms of S from Eq. (6) and then integrates (12) over the 
closed cross section. In view of the fact that V must be a univalued function this leads 


to the following relation: 


s ds oi £ dz (13) 


, Gtr . a 


or, with S from (5), 
k ds dz . . 
= ££ — 5 (14) 
P~ F Gp I 5 
where according to (9) the change in length per winding of a spring has a value 27k. 
Equation (14) is the counterpart of the well-known Bredt formula 6/T = [£(ds/Gt)]/4A? 
for the twist-torque ratio for closed thin-walled sections. 
Examples. We take for a first example the case of a tube with circular cross section 
and with uniform wall thickness. The equations of the center line of the tube wall are 


taken in the form 





r=R-+acosy, z=asiny. (15) 
The integrals occurring in Eqs. (5) and (14) become 
§% _ a [" —_ cospdy 
yr #F. [1 + (a/R) cos y)° 
(16) 
_ _2r (a R)’ 
7 R (1 — (a/R)*}*” 
and 
$ ds a [" : _ dy ae 
Jr RJ, [1 + (a/R) cos y)” 
(17) 
_ma 2+(a R)° 
R* (1 — (a/R)*]°” ° 
The ratio (5) between stress resultant S and applied force P is then 
s > = /PD\ 2713/2 
= jf, Maen (18) 
P = 2xa’ [1 + a/R cos y] 
As the ratio a/R tends to zero the tube stress resultant S approaches the value S = 
PR/2ra” = PR/2A which coincides, as it should, with the value predicted by pure 


torsion theory for an applied couple of magnitude PR. 
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Che maximum stress occurs when cos ¥ = —1 and is of the following magnitude 


Snax _ _R_ [1 — (a/R) 
P Qra’ |i (a/R)|? 
(19) 
. a 3 @) | 
? , 2 sz a al ee é 
2ra’ R t 9\R { 
The deflection-force ratio k/P as given by Eq. (14) becomes 
k_o1 RY [ ~ (2) "|: 4! (2) | 
P~ Gt 2xra* L R 5 \p 
(20) 


i x E 7 g (2) 1 ©) a. | 
Gt 2ra® 8\R 8\R j 


The absence of a term of the form (a/R)’ inside the bracket indicates that the influence 
of the factor in brackets is quite small in all but extreme cases. 

As a second example we consider a tube with rectangular cross section and uniform 
thickness. We designate by R the distance of the center of the cross section from the 
axis, by 2a the width of the tube and by 2b the height of the tube. We find that 


dz 8ab l s 
f Po R*(1 — (a /R)*|? (21) 
and 
ds_ 4(a + b) 1 + (8b — a)(a/k)’ /(b + a) “ 
f ie R? [1 — (a/R)*}’ (22) 


From (21) and (5) the maximum value of the stress resultant S follows in the form 


- R [1 — (a/k)*} 


P ~ 8ab [1 — (a/R)! 


R a a\* 
. mak fis tl (7) a | 


It is noteworthy that the factor in brackets in (23) is independent of the height 2b of 
the tube cross section, and that this factor is somewhat smaller than the corresponding 


factor of Eq. (19) for the circular tube. 
Combination of (14), (21) and (22) gives for the deflection-force ratio of the rec- 





(23) 


tangular tube the following result 
k 1 (a+ dR E 7 (2) [1 4 3b=a (s)'| 
P~ Gt 16a*b’ R b+a\R 


_1@+bR ! 1420>4 G] _ 3b-a (2)'| 
~ Gt 16a°b’ “b+a\R b+a \R/ J 
A comparison with the corresponding formula (20) shows that the square tube shares 


with the circular tube the property that no terms with (a/R)’ occur within the last 
bracket. In contrast to this, terms with (a/R)* do occur whenever one of the sides of 








(24) 
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the rectangle is longer than the other. It may be of interest to list the following special 
cases. Denoting the factor in brackets by 1 + 6, we find that when 


a <b, 1+6~1 +4 2(a/R)*® — 3(a/R)*, (25a) 


and when 


b <a, 1+6~1 — 2(a/R)’ + (a/R). (25b) 


: - ° » /oOF 2 
It may be noted when a < b then 6 assumes a maximum value of 5/27 for (a/R)* = 


1/3, whereas when b < a then 4 is always negative. 


CORRECTIONS* TO THE PAPER 


ON A CLASS OF SINGULAR INTEGRAL EQUATIONS OCCURRING IN PHYSICS 
QUARTERLY OF APPLIED MATHEMATICS 6, 443-448 (1949) 
By H. P. THIELMAN (Jowa State College) 


The limits on the integral in Eq. (B), p. 445 have been omitted. They should have 


been indicated as 0 and ©. 
Equation (a) of Theorem I, p. 445 should read kf(0) — f’(0) = 0 and not kf(0) — 
f"(0) 0 as stated. It should have been stated that f’’(z) in Theorem I, and f‘*(z) in 


Theorem II are assumed to be of order o(e’*) as x goes to infinity. 


*Received June 6, 1949. 
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Proceedings of a symposium on large-scale digital calculating machinery. Jointly sponsored 
by the Navy Department Bureau of Ordnance and Harvard University at the 
Computation Laboratory. Harvard University Press, Cambridge, Massachusetts, 
1948. xxix + 302 pp. $10.00. 


This is a collection of papers and discussions of papers presented at a symposium on large-scale digital 
computing machinery held at Harvard University on January 7-10, 1947. The meeting was sponsored 
jointly by the Navy Department Bureau of Ordnance and Harvard University. The book contains 
numerous photographs and drawings. The technical addresses covered eight sessions dealing with the 
general topics of “Existing Calculating Machines”, “The Logic of Large Scale Calculating Machinery”, 
“Storage Devices”, “Numerical Methods and Suggested Problems for Solution’”’, ‘Sequencing, Coding 
and Problem Preparation’, ‘Input and Output Devices’’, “Conclusions and Open Discussion”’. The state 
of the art seems to have been well surveyed. 

Roun TRUELL 
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Physics in industry: The measurement of stress and strain in solids. Institute of Physies, 
London, 1948. x + 114 pp. $4.00. 


This small book is based on the Proceedings of a Conference held in July 1946 by the Manchester 
and District Branch of the Institute of Physics. The topics covered include wire-resistance and other 
electrical strain gages and instrumentation, acoustic gages, two and three dimensional photoelasticity, 
and X-ray diffraction. Most of the papers are mainly summaries of the highlights of existing information 
and so are of more interest to British than to American workers who have the Proceedings of the Society 
for Experimental Stress Analysis and several books on some phases of the field available. However, there 
are a number of fundamental points concerning the physical behavior of wire-resistance gages which are 
raised and answered by Mr. Eric Jones in his very interesting and important contribution which should 
be read by all. A few of the other features which have not received much attention here are high frequency 
strain gages, the acoustic gage, and the precision setting of distance from film to specimen instead of using 


a reference metal in X-ray back reflection. 
D. C. DrucKER 


Heat conduction with engineering and geological applications. By L. R. Ingersoll, O. J 
Zobel, and A. C. Ingersoll. McGraw-Hill Book Company, Inec., New York and 


London, 1948. xii + 278 pp. $4.00. 


The subject of heat conduction is treated in two ways in various texts. As part of the general subject 
of heat transfer, heat conduction is treated from a practical point of view in various books on engineering 
applications of heat transfer. The mathematical theory, on the other hand, has been treated extensively 
and rigorously in several good texts. The intermediate region which gives sufficient mathematical back- 
ground for a real understanding of the problems involved and yet a sufficient discussion of practical 
applications so as to make obvious the mode of application of the theory has been occupied by the book 
by Ingersoll and Zobel for many years. 

The present work is a thorough revision of the older text with further engineering and geological 
applications. No attempt is made to give a rigorous mathematical foundation for the subject. On the 
other hand, several of the more advanced methods of solution of steady and non-steady heat flow prob- 
lems are presented in a form which makes their ready use possible. Problems of increasing difficulty are 
treated through the book in the following order: Steady-state one-dimension, steady-state more than one 
dimension, periodic heat flow in one dimension, general non-steady heat flow in one dimension, and general 
non-steady heat flow in more than one dimension. A chapter is devoted to the classical treatment of the 
more complex problem of the two conducting media with phase change at the interface as in the forma- 
tion of ice. Two final chapters deal with numerical and graphical methods for solving various types of 
heat conduction problems and a brief discussion of the methods of measuring the various thermal proper- 
ties. Although the major applications are taken from geological problems, the book is a good one regard- 
less of the field of application in which the reader is interested. The authors state as their aims, ‘‘. . . to 
develop the subject with special reference to the needs of the student who has neither time nor mathe- 
matical preparation to pursue the studies at greater length...” and “. . . to point out the many applica- 
tions of which the results are susceptible. . .’’ The book succeeds admirably in living up to these stated 
purposes and is well worth consideration by anyone interested in applications of heat conduction. 


H. W. Eimmons 


Applied mathematics for engineers and scientists. By 8. A. Schelkunoff. D. Van Nostrand 
Co., Inc., Toronto, New York and London, 1948. xi + 472 pp. $6.50. 


This latest addition to the rapidly increasing number of books on ‘‘Applied Mathematics for the 
Engineer, etc.”’, seems particularly appropriate for those scientists whose training has not been that of a 
mathematician. The more elementary material includes a presentation of the notions of differentiation, 





ee 
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integration, and the algebra of complex numbers and vectors. The former topics are not treated at great 
length but no knowledge of the calculus is taken for granted. Chapters are devoted to interpolation, 
solution of algebraic and transcendental equations, power series, and the introduction of exponential and 
related functions. For reference use a chapter on coordinate systems (listing the conventional operators 
in these systems) is included as well as chapters on Bessel functions, Fresnel integrals, and other trans- 
cendental functions. An extensive chapter on second order ordinary differential equations contains 
methods of solution ranging from the most elementary to those for obtaining asymptotic developments. 
The material on partial differential equations stresses the product-series solution representation. Con- 
formal mapping (largely a discussion of certain explicit mappings), the use of the Laplace transform, 
and contour integration are treated, and the book ends with the derivation, from many physical problems, 
of the appropriate differential equations. Nearly all procedures (from integration to approximation 
techniques) are related to physical processes. Those analogs are usually electrical in character and the 
non-electrically trained may well prefer the mathematical method per se. 

Each new topic is introduced in an extremely simple and clear manner at the minor expense, for 
example, of scattering bits of complex variable theory over several chapters. Much of the important 
material in the topies given above has been presented, and the book should be valuable to those interested 
in both introductory material and the conventional techniques applicable to many physical problems. 


G. F. CARRIER 


Computation curves for compressible fluid problems. By. C. L. Dailey and F. C. Wood. 
John Wiley & Sons, Ine., New York and Chapman « Hall, Ltd., London. x + 33 pp. 
$2.00. 

This set of computation curves is intended as a supplement to the text of Liepmann and Puckett 
“Introduction to Aerodynamics of a Compressible Fluid’. Following a brief exposition concerning the 
equations used, there are three sets of charts which are frequently useful in speeding up computational 
work. The first set of graphs are plots of the conventional point functions of the Mach Number (e.g. 
Pressure /stagnation pressure v.s. Mach No., ete.). for example, (1) the shock angle, and (2) the stagna- 
tion pressure ratio, are plotted as functions of incoming Mach number and deflection angle. In section 3 
the corresponding information on conical shocks is given. In the latter section, y is taken as 1.405. In 


the others y = 1.4 
G. F. CARRIER 


Tables of the Bessel functions of the first kind of orders forty through fifty-one. By The Staff 
Of The Computation Laboratory. Harvard University Press, Cambridge, 1948. 620 pp. 
$10.00 
The Bessel Functions Jso(r), Ja(x), «** , Js:(r), are tabulated to ten decimal places. The argument 


varies in steps of .01 forO < x < 100. 
G. F. CARRIER 


Contributions to applied mechanics (Reissner Anniversary Volume). Edited by The Staff 
Of The Department Of Aeronautical Engineering and Applied Mechanics Of The 
Polytechnic Institute of Brooklyn. J. W. Edwards, Ann Arbor, Michigan, 1949. 
viii + 493 pp. $6.50. 

This testimonial to the great contributions of Professor H. Reissner to aeronautical and structural 
engineering, applied mathematics, mechanics, and to physics contains excellent original and expository 
papers by many of his eminent colleagues and friends. Lack of space prevents the inclusion of even the 
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titles of all the papers and permits only the listing of the authors under the general classification of their 
subject. AERODYNAMICS: S. Bergman, R. Paul Harrington and Paul A. Libby, Henry G. Lew, 
D. P. Riabouchinsky, Walter Tollmien. DYNAMICS: Martin Goland, Paul Lieber and M. E. Hamilton, 
Rufus Oldenberger, D. Williams, 8. W. Yuanand M. Morduchow. ELASTICITY AND STRUCTURES: 
L. H. Donnell, K. O. Friedrichs, R. Gran Olsson, Eric Reissner, A. Schleusner, George Schnadel, J. J. 
Stoker, N. J. Hoff, V. L. Salerno, Harold Liebowitz, Bruno A. Boley, Sebastian V. Nardo. ELEC- 
TRICITY: Ronald M. Foster, Reinhold Riidenbergz. MATHEMATICAL METHODS: Hilda Geiringer, 
R. Grammel, Alexander Weinstein. PLASTICITY: R. v. Mises, A. Nadai, Folke K. G. Odqvist. 
PROPULSION: Theodore von Karman, Paul Torda. 
D. C. Drucker 


Tables of generalized sine- and cosine-integral functions: Part I. By The Staff Of The 
Computation Laboratory. Harvard University Press, Cambridge, 1949. vii + 462 
pp. $10.00. 


The twenty-eight introductory pages of this volume are concerned with the definitions of the 
functions to be tabulated, a discussion of the computation methods, the interpolation technique, a brief 
application list with bibliography, and tables of certain coefficients used in computing the tables which 
comprise the remainder of the book. Six functions are computed: They are fo F(a, s)ds where the six 
integrands in question are u! sin u, u“1(1 — cos u), wu! sin u sin 2, u“! sin uw cos 2, uw cos u sin xz, um! 
cos u(1 — cos z), with u(z, a) = (2? + a?)"2, 

These functions are tabulated in the range 0 < a < 1 in steps of .01 forO0 < x < 1 with increments 


s 


Ax = .01. When ais a multiple of .05, they are also given for 1 < x < 2 with Az = 02 and2 <2 <5 


- 


with Ar = .05. For certain other values of a the tables for 1 < x < 2 are also included. For 1 < a < 2, 
the increments are in general larger but not uniform and values of the functions are available for 
0 < 2 < 25 when ais a multiple of .1. 

G. F. CARRIER 


Tables of generalized sine- and cosine-integral functions: Part IT. By The Staff Of The 
Computation Laboratory. Harvard University Press, Cambridge, 1949. 560 pp. $10.00. 


This volume is a continuation of the tables described in the foregoing review. The range of ais 
2<a < 25,0 < 2 < 25. The increments in a are: .05 fora < 5, .1 fori < a < 10, .2 fora > 10. 


~ 


The increments in z are also: .05 for z < 5, .1 for xz < 10, .2 for x < 25. However, the tables run to 
x = 25 only for intermittent values of a. For some values of a they terminate at x = 5, for others at 
z = 10. 

G. F. Carrier 


Vectorial mechanics. By E. A. Milne. Interscience Publishers Inc., New York, 1948. 
xiii + 382 pp. $7.50. 


This text furnishes a thorough and unified treatment of mechanics by Gibbs’ vector and dyad 
analysis. The author’s decision to use only this tool, limits the scope of the text to some extent. Thus, the 
following topics are omitted: (1) the Lagrange generalized coordinates; (2) variational method; (3) 
the integration of the rotational equations of motion of a rigid body in terms of the Eulerian angles, 
which are defined in the text. Further, the author’s very complete introduction to Gibbs’ vector analysis 
leads to the use of such uncommon terms as, “‘tensor of a vector’, ‘‘vector of a tensor’. However, these 
disadvantages are more than compensated for by the author’s clear and unified treatment of topics. 
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llowing discussions should be of interest to the reader: (1) a vector technique for 
solving linear vector differential equations with constant coefficients; (2) an elegant vector treatment of 
Euler’s theorem on rotations about a point in a rigid body; (3) a simple and interesting vector decom- 
he triple vector expansion; (4) a very general tensor formulation of the Gauss and 
s. The author’s treatment of the last mentioned topic is facilitated by introducing 


In particular, th 


position proof of t 
Stokes integral theorem 
the index notation and the basic concepts of tensor algebra. 

The following topics are discussed in the order listed; (1) vector and tensor analysis; (2) systems of 
line vectors with applications to statics and small displacements of a rigid body; (3) kinematics and 
dynamics of a particle and a rigid body. The material of the text is fully illustrated by numerous prob- 
lems. Many of these problems are of intrinsic interest and are completely solved by the author. In the 
general theory and the problems, the author’s methods illustrate the power of the vector approach to 


mechanics. 
In short, the author has written an excellent text, which even the expert may profitably peruse. 


N. Cospurn 


Advances in applied mechanics. Edited by Richard von Mises and Theodore von Karman, 
Academic Press Ine., New York, 1948. viii + 293 pp. $6.80. 


This book is the first of a series intended to present the results of current research in various fields 
of Applied Mechanics in the form of collections of expository monographs and summaries. The parts of 
this volume are described separately. 

Hugh L. Dryden, “Recent advances in the mechanics of boundary layer flow’’. The subject is discussed 
with particular reference to the problem of the stability of a laminar boundary layer and to the nature of 
turbulent flow in boundary layers. The nature of the laminar boundary layer and the effects of boundary 
layer suction are also discussed. An extensive bibliography is presented. 

N. Minorsky, ‘“Modern trends in non-linear mechanics”. After a brief history of the development of 
non-linear mechanics the modern approach is presented. The first section deals with topological methods 
and gives the results obtained and the inherent limitations of this attack. The second section treats the 
analytical methods of Poincaré, van der Pol, Kryloff and Bogoliuboff, and others. The third section deals 
with non-linear resonance and associated phenomena. 

C. B. Biezeno, “Survey of papers on elasticity published in Holland 1940-1946”. As indicated by the 
title this contribution is in the nature of a review. The works discussed are diverse and lie in all of the 
subj ct fields of « lasticity. 

J. M. Burgers, ‘A mathematical model illustrating the theory of turbulence”’. A simplified mathematical 
model is devised with the intention of simulating the problem of turbulence in an incompressible fluid. 
The model does not have the same complicated geometrical character as the turbulence problem but 
retains the essential non-linearity and, in analogue, the characteristics of energy transfer and dissipation. 
From the mathematical treatment results analogous to those of the isotropic theory of turbulence are 
obtained and suggestions are made as to other characteristics of turbulent fluid motion. 

H. Geiringer, ‘(On numerical methods in wave interaction problems’’. The problem of calculating the 
unsteady one-dimensional flow fields in a perfect gas is investigated. The principal treatment is based on 
the Hugoniot shock relations; the approach of J. von Neumann is also presented. 

R. von Mises and M. Schiffer, “On Bergman’s integration method in two-dimensional compressible fluid 
flow’’. Part I presents the general theory of S. Bergman based upon the equation for the stream function 
in the hodograph plane. Part II presents the particular treatment where a simplified pressure-density 
relation is assumed. 

This book and future volumes of the same series make possible the publication of papers of exposi- 
tory, review, or developmental nature in monograph length. This is particularly valuable because of the 
present dearth of monograph journals and the shortage of space in the regular journals in applied me- 
chanics and mathematics. The future of the venture will, of course, depend upon demand and upon the 
quality of profferred contributions. 

To quote from the preface—“The present volume might be considered as indicative of the topics 
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to be dealt with in future issues and as an example for the diversified kinds of approach we wish to culti- 
vate. In both respects, however, the Editors reserve a certain freedom of choice. Suggestions are invited 


and the offering of contributions will be appreciated.” 
Watiace D. Hayes 


Calcolo tensoriale e applicaziont. By Bruno Finzi and Maria Pastori. Nicola Zanichelli 
Editore, Bologna, 1949. vii + 427 pp. Lire 2000. 


The authors’ purpose, as stated in the introduction, is to present the concepts and principal methods 
of the tensor calculus, in order to facilitate its application by mathematicians, physicists, and engineers. 
The reader’s minimum mathematical background is presupposed to be a (first course) knowledge of the 
differential and integral calculus. The exposition is careful and orderly, proceeding by easy stages; and 
the authors have indeed achieved their aim. 

The book consists of an introduction, ten chapters, and a bibliography. The introduction is in the 
nature of an orientation; the reader is “‘briefed’’ on the ground to be covered and acquainted with the 
concepts that will occur again and again, e.g., vectors, tensors, contravariant superscripts, covariant 
subscripts, etc. A pleasant feature of the book is the table of principal formulas (‘‘formule notevoli’’) 
which follows every chapter. The first three chapters deal with the tensor calculus proper: vector fields, 
tensors and tensor algebra, and “‘omografie vettoriali’”’ (linear vector functions). The next four chapters 
are concerned with applications to differential geometry: tensor fields in Euclidean and non-Euclidean 
spaces, surface geometry and Riemannian geometry. The last three chapters contain the applications to 
mathematical physics referred to in the first paragraph of the introduction. The authors have succeeded 
in compressing an unusually large amount of material into comparatively few pages. The chapter headings 
are: VIII Mechanics of Deformable Continua, IX Electromagnetic Theory, and X Relativity Theory. 
The book concludes with a bibliography of representative works on various phases of the subject. 


J. B. Diaz 


Supersonic flow and shock waves. By R. Courant and K. O. Friedrichs. Interscience Pub- 
lishers, Inc., New York and London, 1948. xvi + 464 pp. $7.00. 


This book is primarily concerned with the hyperbolic problems associated with the behavior of con- 
tinua. Except for brief remarks about shallow water waves and discontinuous waves in plastic media, 
attention is confined to gas flows. The opening chapter provides the fundamental notions (thermody- 
namic concepts, and basic mechanical laws) which form the basis of the analysis of these phenomena. 
This is followed by a chapter on the mathematics of hyperbolic systems of second order. Unsteady one- 
dimensional problems are then considered; shock interaction, detonation, and plastic waves are dis- 
cussed in this section. Chapter IV deals with steady plane flow; it utilizes the hodograph method, to 
present certain isentropic flows, and contains the analysis of shock interactions, and treats the flow past 
obstacles, ete. by the perturbation method. The book ends with briefer discussions of nozzle flows, jets, 
axially symmetric flows, spherical waves and certain conical flow problems. 

The book seems well adapted for the presentation to students of that portion of compressible fluid 
theory which is essentially hyperbolic. It is obviously not intended as a self contained text on the general 


field of compressible flow theory. 
G. F, CARRIER 























